Engineering Logical Inflammation Sensing Circuit for Gut Modulation

BioCRNPyler Models and Code for Fig 5C, 6BCD

Liana Merk

In [1]:
from bioscrape.types import Model
from bioscrape.simulator import py_simulate_model

from biocrnpyler import *

import numpy as np
import pylab as plt
import tqdm
import pandas as pd

import matplotlib as mpl
import holoviews as hv
import bokeh.plotting
import bokeh.io
from bokeh.models import Range1d
from bokeh.layouts import gridplot
bokeh.io.output_notebook()
hv.extension('bokeh')

from bokeh.models import PrintfTickFormatter

import seaborn as sns
sns.set()
%matplotlib inline
Loading BokehJS ...

Tetrathionate Two Component System

In [2]:
#Parameters": List of tuples [("paramter name" (string)": value (float))]
parameters = {"k1b": 100, "k1u": 1, 
              "k2b": 0.3, "k2u": 10, 
              "k3b": 0.3, "k3u": 10, 
              "k4b": 1.6, "k4u": 0.016,
              "k5b": 0.0001, "k5u":1, 
              "k6b": 0.0001, "k6u": 0.1, 
              "k7b": 1, "k7u": 1, 
              "k8b": 0.0083, "k8u": 5, 
              "k9b": 0.3, "k9u": 3,
              "k10b": 0.1, "k10u":1, 
              "k11b": 100, "k11u": 1,     
              "ktx": 0.09, "ktl": 0.05,
              "k_dephos": 50, "delta": .9}

Define all species and reactions

Some Basic Species

In [3]:
# Ribosome
R = Species("R")

# RNA Polymerase
P = Species("P")

Reaction 1: Polymerase binding to P1d

In [4]:
# Consitutive Promoter P1d
P1d = Species("P1d")

# Polymerase bound to P1d
P_P1d = Species("P:P1d")

# P + P1d <-> P1d:P
R1 = Reaction([P, P1d], [P_P1d], k = parameters['k1b'], k_rev = parameters['k1u'])

Transcription 1 (ttrS and ttrR)

In [5]:
# ttrS Transcript
ttrS_T = Species("ttrS_T")

# ttrR Transcript
ttrR_T = Species("ttrR_T")

# P1d:P -> ttrS_T + ttrR_T
Rtx1 = Reaction([P_P1d], [ttrS_T, ttrR_T, P, P1d], k = parameters['ktx'])

Reaction 2 and translation 1: ttrS

In [6]:
# ttrS trasncript bound to Ribosome
ttrS_T_R = Species("ttrS_T:R")

# ttrS
ttrS = Species("ttrS")

# ttrS_T + R <-> ttrS_T:R
R2 = Reaction([ttrS_T, R], [ttrS_T_R], k = parameters['k2b'], k_rev = parameters['k2u'])

# ttrS Translation
Rtl1 = Reaction([ttrS_T_R], [ttrS_T, ttrS, R], k = parameters['ktl'])

Reaction 3 and translation 2: ttrR

In [7]:
# ttrR trasncript bound to Ribosome
ttrR_T_R = Species("ttrR_T:R")

# ttrR
ttrR = Species("ttrR")

# ttrR_T + R <-> ttrR_T:R
R3 = Reaction([ttrR_T, R], [ttrR_T_R], k = parameters['k3b'], k_rev = parameters['k3u'])

# ttrR Translation
Rtl2 = Reaction([ttrR_T_R], [ttrR_T, ttrR, R], k = parameters['ktl'])

Reaction 4: Tetrathionate binding to ttrS, Causing Phosphorylation

In [8]:
# Tetrathionate
tt = Species("tt")

# Phosphorylated ttrS
ttrS_P = Species("ttrS_P")

# ttrS + tt <-> ttrS_P + tt
R4 = Reaction([ttrS, tt], [ttrS_P, tt], k = parameters['k4b'], k_rev = parameters['k4u'])

Reaction 5: ttrR binding to dephosphorylated ttrS

In [9]:
# ttrR bound to ttrS
ttrR_ttrS = Species("ttrR_ttrS")

# ttrR + ttrS <-> ttrR:ttrS
R5 = Reaction([ttrS, ttrR], [ttrR_ttrS], k = parameters['k5b'], k_rev = parameters['k5u'])

Reaction 6: ttrR binding to phosphorylated ttrS

In [10]:
# ttrR bound to phosphorylated ttrS
ttrR_ttrS_P = Species("ttrR_ttrS_P")

# ttrR + ttrS <-> ttrR:ttrS:P
R6 = Reaction([ttrS_P, ttrR], [ttrR_ttrS_P], k = parameters['k6b'], k_rev = parameters['k6u'])

Reaction 7: Phosphorylation transfer ttrS to ttrR

In [11]:
# Phosphorylated ttrR
ttrR_P = Species("ttrR_P")

# ttrR:ttrS:P <-> ttrR:P + ttrS
R7 = Reaction([ttrR_ttrS_P], [ttrR_P, ttrS], k = parameters['k7b'], k_rev = parameters['k7u'])

Dephosphorylated

In [12]:
# ttrR:P -> ttrR
rxn_dephos = Reaction([ttrR_P], [ttrR], k = parameters['k_dephos'])

Reaction 8: Dimerization of phosphorylated ttrR

In [13]:
# ttrR:P_dimer
ttrR_P_dimer = Species("ttrR:P_dimer")

# 2ttrR:P <-> ttrR:P_dimer
R8 = Reaction([ttrR_P, ttrR_P], [ttrR_P_dimer], k = parameters['k8b'], k_rev = parameters['k8u'])

Reaction 9: Phosphorylated ttrR dimer activated promoter pttr0, turns to pttr1

In [14]:
# Inactive pttr
pttr0 = Species("pttr0")

# Inactive pttr
pttr1 = Species("pttr1")

# ttrR:P_dimer + pttr0 <-> pttr1
R9 = Reaction([ttrR_P_dimer, pttr0], [pttr1], k = parameters['k9b'], k_rev = parameters['k9u'])

Reaction 10: Polymerase binding to active promoter

In [15]:
# Active pttr to polymerase
pttr1_P = Species("pttr1_P")

# P + pttr1 <-> pttr1:P
Rtx2 = Reaction([P, pttr1], [pttr1_P], k = parameters['k10b'], k_rev = parameters['k10u'])

Transcription 2: active pttr produces GFP transcript

In [16]:
# Active pttr to polymerase
GFP_T = Species("GFP_T")

# P + pttr1 <-> pttr1:P
R10 = Reaction([pttr1_P], [GFP_T, P, pttr1], k = parameters['ktx'])

Reaction 11: Ribosome binding to GFP transcript and translation 3

In [17]:
# Active pttr to polymerase
GFP_T_R = Species("GFP_T_R")

# GFP
GFP = Species("GFP")

# P + pttr1 <-> pttr1:P
R11 = Reaction([GFP_T, R], [GFP_T_R], k = parameters['k11b'], k_rev = parameters['k11u'])

# ttrS Translation
Rtl3 = Reaction([GFP_T_R], [GFP, GFP_T, R], k = parameters['ktl'])

Degradation Reactions

In [18]:
# RNA deg
rxnd_ttrS_T = Reaction([ttrS_T], [], k = parameters['delta'])
rxnd_ttrR_T = Reaction([ttrR_T], [], k = parameters['delta'])
rxnd_GFP_T = Reaction([GFP_T], [], k = parameters['delta'])

# Protein Deg
rxnd_ttrS = Reaction([ttrS], [], k = parameters['delta'])
rxnd_ttrR = Reaction([ttrR], [], k = parameters['delta'])
rxnd_GFP = Reaction([GFP], [], k = parameters['delta'])

Now, pull all the reactions together

In [19]:
reactions = [R1, R2, R3, R4, R5, R6, R7, R8, R9, R10, R11,
             Rtx1, Rtx2,
             Rtl1, Rtl2, Rtl3,
             rxnd_ttrS_T, rxnd_ttrR_T, rxnd_GFP_T,
             rxnd_ttrS, rxnd_ttrR, rxnd_GFP]
In [20]:
species = [R, P, P1d, P_P1d, ttrS_T, ttrR_T, ttrR_T_R, ttrS_T_R,
           tt, ttrS_P, ttrR_ttrS, ttrR_ttrS_P, ttrR_P, ttrR_P_dimer,
           pttr0, pttr1, pttr1_P, GFP_T, GFP_T_R, GFP]
In [21]:
# Create CRN
CRN_tt = ChemicalReactionNetwork(species = species, reactions = reactions)

print(CRN_tt.pretty_print(show_rates = True))
Species (20) = {0. R, 1. P, 2. P1d, 3. P:P1d, 4. ttrS_T, 5. ttrR_T, 6. ttrR_T:R, 7. ttrS_T:R, 8. tt, 9. ttrS_P, 10. ttrR_ttrS, 11. ttrR_ttrS_P, 12. ttrR_P, 13. ttrR:P_dimer, 14. pttr0, 15. pttr1, 16. pttr1_P, 17. GFP_T, 18. GFP_T_R, 19. GFP}
Reactions (22) = [
0. P + P1d <--> P:P1d        
        massaction: k_f(P,P1d)=100*P*P1d
        k_r(P:P1d)=1*P:P1d
1. ttrS_T + R <--> ttrS_T:R        
        massaction: k_f(ttrS_T,R)=0.3*ttrS_T*R
        k_r(ttrS_T:R)=10*ttrS_T:R
2. ttrR_T + R <--> ttrR_T:R        
        massaction: k_f(ttrR_T,R)=0.3*ttrR_T*R
        k_r(ttrR_T:R)=10*ttrR_T:R
3. ttrS + tt <--> ttrS_P + tt        
        massaction: k_f(ttrS,tt)=1.6*ttrS*tt
        k_r(ttrS_P,tt)=0.016*ttrS_P*tt
4. ttrS + ttrR <--> ttrR_ttrS        
        massaction: k_f(ttrS,ttrR)=0.0001*ttrS*ttrR
        k_r(ttrR_ttrS)=1*ttrR_ttrS
5. ttrS_P + ttrR <--> ttrR_ttrS_P        
        massaction: k_f(ttrS_P,ttrR)=0.0001*ttrS_P*ttrR
        k_r(ttrR_ttrS_P)=0.1*ttrR_ttrS_P
6. ttrR_ttrS_P <--> ttrR_P + ttrS        
        massaction: k_f(ttrR_ttrS_P)=1*ttrR_ttrS_P
        k_r(ttrR_P,ttrS)=1*ttrR_P*ttrS
7. 2 ttrR_P <--> ttrR:P_dimer        
        massaction: k_f(ttrR_P)=0.0083*ttrR_P^2
        k_r(ttrR:P_dimer)=5*ttrR:P_dimer
8. ttrR:P_dimer + pttr0 <--> pttr1        
        massaction: k_f(ttrR:P_dimer,pttr0)=0.3*ttrR:P_dimer*pttr0
        k_r(pttr1)=3*pttr1
9. pttr1_P --> GFP_T + P + pttr1        
        massaction: k_f(pttr1_P)=0.09*pttr1_P
10. GFP_T + R <--> GFP_T_R        
        massaction: k_f(GFP_T,R)=100*GFP_T*R
        k_r(GFP_T_R)=1*GFP_T_R
11. P:P1d --> ttrS_T + ttrR_T + P + P1d        
        massaction: k_f(P:P1d)=0.09*P:P1d
12. P + pttr1 <--> pttr1_P        
        massaction: k_f(P,pttr1)=0.1*P*pttr1
        k_r(pttr1_P)=1*pttr1_P
13. ttrS_T:R --> ttrS_T + ttrS + R        
        massaction: k_f(ttrS_T:R)=0.05*ttrS_T:R
14. ttrR_T:R --> ttrR_T + ttrR + R        
        massaction: k_f(ttrR_T:R)=0.05*ttrR_T:R
15. GFP_T_R --> GFP + GFP_T + R        
        massaction: k_f(GFP_T_R)=0.05*GFP_T_R
16. ttrS_T -->         
        massaction: k_f(ttrS_T)=0.9*ttrS_T
17. ttrR_T -->         
        massaction: k_f(ttrR_T)=0.9*ttrR_T
18. GFP_T -->         
        massaction: k_f(GFP_T)=0.9*GFP_T
19. ttrS -->         
        massaction: k_f(ttrS)=0.9*ttrS
20. ttrR -->         
        massaction: k_f(ttrR)=0.9*ttrR
21. GFP -->         
        massaction: k_f(GFP)=0.9*GFP
]
In [22]:
# An initial condition for each species (uninitialized species default to 0)
x0 = {
    "tt":200,
    "P1d":5,
    "pttr0":5,
    "P":25,
    "R":100,
}

timepoints = np.linspace(0, 600, 100)
In [23]:
Results = CRN_tt.simulate_with_bioscrape(timepoints,x0)
/Users/lianamerk/anaconda3/envs/python3.7-bioscrape/lib/python3.7/site-packages/biocrnpyler-0.1-py3.7.egg/biocrnpyler/chemical_reaction_network.py:870: UserWarning: The following species are uninitialized and their value has been defaulted to 0: ttrS, ttrR, 
In [24]:
Results.head()
Out[24]:
R P P1d P:P1d ttrS_T ttrR_T ttrR_T:R ttrS_T:R tt ttrS_P ... ttrR:P_dimer pttr0 pttr1 pttr1_P GFP_T GFP_T_R GFP ttrS ttrR time
0 100.000000 25.000000 5.000000 0.000000 0.000000 0.000000 0.000000 0.000000 200.0 0.000000 ... 0.000000e+00 5.0 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000 0.000000 0.000000
1 97.824944 20.002723 0.002723 4.997277 0.375266 0.375266 1.087528 1.087528 200.0 0.194098 ... 1.163234e-15 5.0 3.084548e-16 2.295029e-16 -1.763695e-18 -1.645280e-14 -9.145889e-16 0.002103 0.053622 6.060606
2 97.280408 20.002723 0.002723 4.997277 0.468971 0.468971 1.359796 1.359796 200.0 0.549631 ... 2.082295e-13 5.0 7.285205e-14 8.566282e-14 1.728165e-19 -5.795336e-15 -5.961676e-16 0.005691 0.073832 12.121212
3 97.146497 20.002723 0.002723 4.997277 0.492174 0.492174 1.426751 1.426751 200.0 0.928944 ... 2.437093e-12 5.0 9.621154e-13 1.344522e-12 4.207668e-17 2.741544e-13 9.917796e-15 0.009484 0.078838 18.181818
4 97.113649 20.002723 0.002723 4.997277 0.497875 0.497875 1.443175 1.443175 200.0 1.298786 ... 1.128285e-11 5.0 4.751947e-12 7.217950e-12 3.139122e-16 2.285017e-12 9.528891e-14 0.013174 0.080064 24.242424

5 rows × 23 columns

Holoviews plotting of Bioscrape results

Sort out Plotting

In [25]:
def plot_holoviews_ttr_bioscrape(results_melted):
    plot_group = []

    for row in results_melted.iterrows():
        sp = row[1]['species']

        if sp in ["P1d", "P:P1d", "pttr0", "pttr1", "pttr1_P"]:
            plot_group.append('promoters')

        elif sp in ["ttrS_T", "ttrR_T", "ttrR_T:R", "ttrS_T:R", "GFP_T", "GFP_T_R"]:
            plot_group.append('mRNA')

        elif sp in ["ttrS_P", "ttrR_ttrS", "ttrR_ttrS_P", "ttrR_P", "ttrR:P_dimer"]:
            plot_group.append('complexes')

        elif sp in ["ttrR", "ttrS", "GFP"]:
            plot_group.append('proteins')

        elif sp in ["P",   "R", "tt"]:
            plot_group.append('machinery and signals')
        else:
            print(sp)

    results_melted['plot_group'] = plot_group
    

    plot_list = []
    for group in results_melted['plot_group'].unique():
        plot_list.append(hv.Curve(
            data=results_melted.loc[results_melted['plot_group'] == group],
            kdims=['time', 'value'],
            vdims = ['species']
        ).groupby(
            [ 'species']
        ).opts(
            tools=['hover'],
#             legend_position = 'right',
            show_grid=True,
            title = f'{group}',
            width = 500,
            #shared_axes=False
        ).overlay('species').opts(legend_position = 'right',))
        
    return plot_list

Now, proceed with plotting

In [26]:
results_melted = pd.melt(Results, var_name='species', id_vars=['time'])
results_melted.head()
Out[26]:
time species value
0 0.000000 R 100.000000
1 6.060606 R 97.824944
2 12.121212 R 97.280408
3 18.181818 R 97.146497
4 24.242424 R 97.113649
In [27]:
plot_list = plot_holoviews_ttr_bioscrape(results_melted)
hv.Layout(plot_list).cols(2).opts(shared_axes=False)
Out[27]:

Tuning Function

In [28]:
def create_reactions(parameters):
    # P + P1d <-> P1d:P
    R1 = Reaction([P, P1d], [P_P1d], k = parameters['k1b'], k_rev = parameters['k1u'])

    # P1d:P -> ttrS_T + ttrR_T
    Rtx1 = Reaction([P_P1d], [ttrS_T, ttrR_T, P, P1d], k = parameters['ktx'])
    
    # ttrS_T + R <-> ttrS_T:R
    R2 = Reaction([ttrS_T, R], [ttrS_T_R], k = parameters['k2b'], k_rev = parameters['k2u'])

    # ttrS Translation
    Rtl1 = Reaction([ttrS_T_R], [ttrS_T, ttrS, R], k = parameters['ktl'])
    
    # ttrR_T + R <-> ttrR_T:R
    R3 = Reaction([ttrR_T, R], [ttrR_T_R], k = parameters['k3b'], k_rev = parameters['k3u'])

    # ttrR Translation
    Rtl2 = Reaction([ttrR_T_R], [ttrR_T, ttrR, R], k = parameters['ktl'])
    
    # ttrS + tt <-> ttrS_P + tt
    R4 = Reaction([ttrS, tt], [ttrS_P, tt], k = parameters['k4b'], k_rev = parameters['k4u'])
    
    # ttrR + ttrS <-> ttrR:ttrS
    R5 = Reaction([ttrS, ttrR], [ttrR_ttrS], k = parameters['k5b'], k_rev = parameters['k5u'])

    # ttrR + ttrS <-> ttrR:ttrS:P
    R6 = Reaction([ttrS_P, ttrR], [ttrR_ttrS_P], k = parameters['k6b'], k_rev = parameters['k6u'])

    # ttrR:ttrS:P <-> ttrR:P + ttrS
    R7 = Reaction([ttrR_ttrS_P], [ttrR_P, ttrS], k = parameters['k7b'], k_rev = parameters['k7u'])
    
    # ttrR:P -> ttrR
    rxn_dephos = Reaction([ttrR_P], [ttrR], k = parameters['k_dephos'])

    # 2ttrR:P <-> ttrR:P_dimer
    R8 = Reaction([ttrR_P, ttrR_P], [ttrR_P_dimer], k = parameters['k8b'], k_rev = parameters['k8u'])
    
    # ttrR:P_dimer + pttr0 <-> pttr1
    R9 = Reaction([ttrR_P_dimer, pttr0], [pttr1], k = parameters['k9b'], k_rev = parameters['k9u'])

    # P + pttr1 <-> pttr1:P
    Rtx2 = Reaction([P, pttr1], [pttr1_P], k = parameters['k10b'], k_rev = parameters['k10u'])

    # P + pttr1 <-> pttr1:P
    R10 = Reaction([pttr1_P], [GFP_T, P, pttr1], k = parameters['ktx'])

    # P + pttr1 <-> pttr1:P
    R11 = Reaction([GFP_T, R], [GFP_T_R], k = parameters['k11b'], k_rev = parameters['k11u'])

    # ttrS Translation
    Rtl3 = Reaction([GFP_T_R], [GFP, GFP_T, R], k = parameters['ktl'])

    # RNA deg
    rxnd_ttrS_T = Reaction([ttrS_T], [], k = parameters['delta'])
    rxnd_ttrR_T = Reaction([ttrR_T], [], k = parameters['delta'])
    rxnd_GFP_T = Reaction([GFP_T], [], k = parameters['delta'])

    # Protein Deg
    rxnd_ttrS = Reaction([ttrS], [], k = parameters['delta'])
    rxnd_ttrR = Reaction([ttrR], [], k = parameters['delta'])
    rxnd_GFP = Reaction([GFP], [], k = parameters['delta'])
    
    reactions = [R1, R4, R5, R6, R7, R8, R9, R10, R11,
                 Rtx1, Rtx2,
                 Rtl1, Rtl2, Rtl3,
                 rxnd_ttrS_T, rxnd_ttrR_T, rxnd_GFP_T,
                 rxnd_ttrS, rxnd_ttrR, rxnd_GFP]
    
    return reactions

RBS Tuning in Tetrathionate Circuit

We will change k2u, and k3u, the reaction rates that define ribosome unbinding from ttrS and ttrR transcripts, respectively. We will save the simulations as csv's to plot in another notebook such that the styles match well.

In [29]:
all_sims = pd.DataFrame()
In [30]:
#Parameters": List of tuples [("paramter name" (string)": value (float))]
parameters = {"k1b": 100, "k1u": 1, 
              "k2b": 3, "k2u": 10, 
              "k3b": 0.00003, "k3u": 10, 
              "k4b": 5, "k4u": 0.5,
              "k5b": 5, "k5u":6, 
              "k6b": 5, "k6u": 1, 
              "k7b": 0.01, "k7u": 1, 
              "k8b": 80, "k8u": 5, 
              "k9b": 0.004, "k9u": 3,
              "k10b": 50, "k10u":1, 
              "k11b": 10, "k11u": 1,     
              "ktx": 0.1, "ktl": 1,
              "k_dephos": 7, "delta": .9}
In [31]:
def plot_ttr_RBS_tuning(parameters, x0, k2u, k3u, title, timepoints):
    experimental_ttr_concs = np.array([ 0, 0.1, 0.25, 0.5, 0.75, 1, 5, 10])

    # Get color palette in order
    r = list(bokeh.palettes.viridis(10))
    r.reverse()
    r = r[2:]

    # ttrS_T + R <-> ttrS_T:R
    R2 = Reaction([ttrS_T, R], [ttrS_T_R], k = parameters['k2b'], k_rev = k2u)

    # ttrR_T + R <-> ttrR_T:R
    R3 = Reaction([ttrR_T, R], [ttrR_T_R], k = parameters['k3b'], k_rev = k3u)


    reactions = create_reactions(parameters)
    reactions.append(R2)
    reactions.append(R3)
    
    species = [R, P, P1d, P_P1d, ttrS_T, ttrR_T, ttrR_T_R, ttrS_T_R,
           tt, ttrS_P, ttrR_ttrS, ttrR_ttrS_P, ttrR_P, ttrR_P_dimer,
           pttr0, pttr1, pttr1_P, GFP_T, GFP_T_R, GFP]
    
    # Recompile the CRN with the right parameters
    CRN_curr = ChemicalReactionNetwork(species = species, reactions = reactions)
    
    p = bokeh.plotting.figure(width = 300, height= 250, x_axis_label='Time (hr)', y_axis_label = '[GFP]', title = title)
    i = 0

    for tt_conc in experimental_ttr_concs:
        x0["tt"] = tt_conc
        
        Results = CRN_curr.simulate_with_bioscrape(timepoints,x0)
        all_sims[f'{title}_{tt_conc}'] = Results['GFP']

        p.line(timepoints/60, Results['GFP'],color = r[i], legend_label = f'{tt_conc}', line_width = 2)
        i+=1

    p.legend.location = 'top_left'
    if k2u != 20:
        p.legend.visible = False
        
    p.y_range = Range1d(-0.1,2.5)

    return p
In [32]:
timepoints = np.linspace(0, 1400, 1400)
title_list = ['Very Low (LM15) Simulation', 'Low (LM14) Simulation', 'Medium (LM18) Simulation', 'High (LM19) Simulation']
In [33]:
x0 = {
    "tt":200,
    # PSC101 has copy number 5
    "P1d":5,
    "pttr0":5,
    # Polymerase: 15 uM assuming ~10000 Polymerase molecules / E. coli with a volume 1 um^3
    "P":15,
    # Ribosome: 120 uM assuming ~72000 Ribosomes / E. coli with a volume 1 um^3
    "R":120,
}
In [34]:
parameters = {"k1b": 10, "k1u": 0.0001, 
              "k2b": 0.3, "k2u": 10, 
              "k3b": 0.3, "k3u": 10, 
              "k4b": 1.6, "k4u": 0.016,
              "k5b": 0.0001, "k5u":1, 
              "k6b": 0.0001, "k6u": 0.1, 
              "k7b": 1, "k7u": 1, 
              "k8b": 0.0083, "k8u": 0.5, 
              "k9b": 0.3, "k9u": 3,
              "k10b": 10, "k10u":0.0001, 
              "k11b": 100, "k11u": 1,     
              "ktx": 0.1, "ktl":0.33,
              "k_dephos": 50, "delta": .9}
In [35]:
# Define higher off rates
k2u = [20, 12, 11, 10]
k3u = [15, 8, 6, 5]

plot_list = []

for i in range(len(k3u)):
    plot_list.append(plot_ttr_RBS_tuning(parameters, x0, k2u[i], k3u[i], title_list[i], timepoints))

bokeh.io.show(gridplot([plot_list]))
/Users/lianamerk/anaconda3/envs/python3.7-bioscrape/lib/python3.7/site-packages/biocrnpyler-0.1-py3.7.egg/biocrnpyler/chemical_reaction_network.py:870: UserWarning: The following species are uninitialized and their value has been defaulted to 0: ttrS, ttrR, 

Tetrathionate Two Component System

In [2]:
#Parameters": List of tuples [("paramter name" (string)": value (float))]
parameters = {"k1b": 100, "k1u": 1, 
              "k2b": 0.3, "k2u": 10, 
              "k3b": 0.3, "k3u": 10, 
              "k4b": 1.6, "k4u": 0.016,
              "k5b": 0.0001, "k5u":1, 
              "k6b": 0.0001, "k6u": 0.1, 
              "k7b": 1, "k7u": 1, 
              "k8b": 0.0083, "k8u": 5, 
              "k9b": 0.3, "k9u": 3,
              "k10b": 0.1, "k10u":1, 
              "k11b": 100, "k11u": 1,     
              "ktx": 0.09, "ktl": 0.05,
              "k_dephos": 50, "delta": .9}

Define all species and reactions

Some Basic Species

In [3]:
# Ribosome
R = Species("R")

# RNA Polymerase
P = Species("P")

Reaction 1: Polymerase binding to P1d

In [4]:
# Consitutive Promoter P1d
P1d = Species("P1d")

# Polymerase bound to P1d
P_P1d = Species("P:P1d")

# P + P1d <-> P1d:P
R1 = Reaction([P, P1d], [P_P1d], k = parameters['k1b'], k_rev = parameters['k1u'])

Transcription 1 (ttrS and ttrR)

In [5]:
# ttrS Transcript
ttrS_T = Species("ttrS_T")

# ttrR Transcript
ttrR_T = Species("ttrR_T")

# P1d:P -> ttrS_T + ttrR_T
Rtx1 = Reaction([P_P1d], [ttrS_T, ttrR_T, P, P1d], k = parameters['ktx'])

Reaction 2 and translation 1: ttrS

In [6]:
# ttrS trasncript bound to Ribosome
ttrS_T_R = Species("ttrS_T:R")

# ttrS
ttrS = Species("ttrS")

# ttrS_T + R <-> ttrS_T:R
R2 = Reaction([ttrS_T, R], [ttrS_T_R], k = parameters['k2b'], k_rev = parameters['k2u'])

# ttrS Translation
Rtl1 = Reaction([ttrS_T_R], [ttrS_T, ttrS, R], k = parameters['ktl'])

Reaction 3 and translation 2: ttrR

In [7]:
# ttrR trasncript bound to Ribosome
ttrR_T_R = Species("ttrR_T:R")

# ttrR
ttrR = Species("ttrR")

# ttrR_T + R <-> ttrR_T:R
R3 = Reaction([ttrR_T, R], [ttrR_T_R], k = parameters['k3b'], k_rev = parameters['k3u'])

# ttrR Translation
Rtl2 = Reaction([ttrR_T_R], [ttrR_T, ttrR, R], k = parameters['ktl'])

Reaction 4: Tetrathionate binding to ttrS, Causing Phosphorylation

In [8]:
# Tetrathionate
tt = Species("tt")

# Phosphorylated ttrS
ttrS_P = Species("ttrS_P")

# ttrS + tt <-> ttrS_P + tt
R4 = Reaction([ttrS, tt], [ttrS_P, tt], k = parameters['k4b'], k_rev = parameters['k4u'])

Reaction 5: ttrR binding to dephosphorylated ttrS

In [9]:
# ttrR bound to ttrS
ttrR_ttrS = Species("ttrR_ttrS")

# ttrR + ttrS <-> ttrR:ttrS
R5 = Reaction([ttrS, ttrR], [ttrR_ttrS], k = parameters['k5b'], k_rev = parameters['k5u'])

Reaction 6: ttrR binding to phosphorylated ttrS

In [10]:
# ttrR bound to phosphorylated ttrS
ttrR_ttrS_P = Species("ttrR_ttrS_P")

# ttrR + ttrS <-> ttrR:ttrS:P
R6 = Reaction([ttrS_P, ttrR], [ttrR_ttrS_P], k = parameters['k6b'], k_rev = parameters['k6u'])

Reaction 7: Phosphorylation transfer ttrS to ttrR

In [11]:
# Phosphorylated ttrR
ttrR_P = Species("ttrR_P")

# ttrR:ttrS:P <-> ttrR:P + ttrS
R7 = Reaction([ttrR_ttrS_P], [ttrR_P, ttrS], k = parameters['k7b'], k_rev = parameters['k7u'])

Dephosphorylated

In [12]:
# ttrR:P -> ttrR
rxn_dephos = Reaction([ttrR_P], [ttrR], k = parameters['k_dephos'])

Reaction 8: Dimerization of phosphorylated ttrR

In [13]:
# ttrR:P_dimer
ttrR_P_dimer = Species("ttrR:P_dimer")

# 2ttrR:P <-> ttrR:P_dimer
R8 = Reaction([ttrR_P, ttrR_P], [ttrR_P_dimer], k = parameters['k8b'], k_rev = parameters['k8u'])

Reaction 9: Phosphorylated ttrR dimer activated promoter pttr0, turns to pttr1

In [14]:
# Inactive pttr
pttr0 = Species("pttr0")

# Inactive pttr
pttr1 = Species("pttr1")

# ttrR:P_dimer + pttr0 <-> pttr1
R9 = Reaction([ttrR_P_dimer, pttr0], [pttr1], k = parameters['k9b'], k_rev = parameters['k9u'])

Reaction 10: Polymerase binding to active promoter

In [15]:
# Active pttr to polymerase
pttr1_P = Species("pttr1_P")

# P + pttr1 <-> pttr1:P
Rtx2 = Reaction([P, pttr1], [pttr1_P], k = parameters['k10b'], k_rev = parameters['k10u'])

Transcription 2: active pttr produces GFP transcript

In [16]:
# Active pttr to polymerase
GFP_T = Species("GFP_T")

# P + pttr1 <-> pttr1:P
R10 = Reaction([pttr1_P], [GFP_T, P, pttr1], k = parameters['ktx'])

Reaction 11: Ribosome binding to GFP transcript and translation 3

In [17]:
# Active pttr to polymerase
GFP_T_R = Species("GFP_T_R")

# GFP
GFP = Species("GFP")

# P + pttr1 <-> pttr1:P
R11 = Reaction([GFP_T, R], [GFP_T_R], k = parameters['k11b'], k_rev = parameters['k11u'])

# ttrS Translation
Rtl3 = Reaction([GFP_T_R], [GFP, GFP_T, R], k = parameters['ktl'])

Degradation Reactions

In [18]:
# RNA deg
rxnd_ttrS_T = Reaction([ttrS_T], [], k = parameters['delta'])
rxnd_ttrR_T = Reaction([ttrR_T], [], k = parameters['delta'])
rxnd_GFP_T = Reaction([GFP_T], [], k = parameters['delta'])

# Protein Deg
rxnd_ttrS = Reaction([ttrS], [], k = parameters['delta'])
rxnd_ttrR = Reaction([ttrR], [], k = parameters['delta'])
rxnd_GFP = Reaction([GFP], [], k = parameters['delta'])

Now, pull all the reactions together

In [19]:
reactions = [R1, R2, R3, R4, R5, R6, R7, R8, R9, R10, R11,
             Rtx1, Rtx2,
             Rtl1, Rtl2, Rtl3,
             rxnd_ttrS_T, rxnd_ttrR_T, rxnd_GFP_T,
             rxnd_ttrS, rxnd_ttrR, rxnd_GFP]
In [20]:
species = [R, P, P1d, P_P1d, ttrS_T, ttrR_T, ttrR_T_R, ttrS_T_R,
           tt, ttrS_P, ttrR_ttrS, ttrR_ttrS_P, ttrR_P, ttrR_P_dimer,
           pttr0, pttr1, pttr1_P, GFP_T, GFP_T_R, GFP]
In [21]:
# Create CRN
CRN_tt = ChemicalReactionNetwork(species = species, reactions = reactions)

print(CRN_tt.pretty_print(show_rates = True))
Species (20) = {0. R, 1. P, 2. P1d, 3. P:P1d, 4. ttrS_T, 5. ttrR_T, 6. ttrR_T:R, 7. ttrS_T:R, 8. tt, 9. ttrS_P, 10. ttrR_ttrS, 11. ttrR_ttrS_P, 12. ttrR_P, 13. ttrR:P_dimer, 14. pttr0, 15. pttr1, 16. pttr1_P, 17. GFP_T, 18. GFP_T_R, 19. GFP}
Reactions (22) = [
0. P + P1d <--> P:P1d        
        massaction: k_f(P,P1d)=100*P*P1d
        k_r(P:P1d)=1*P:P1d
1. ttrS_T + R <--> ttrS_T:R        
        massaction: k_f(ttrS_T,R)=0.3*ttrS_T*R
        k_r(ttrS_T:R)=10*ttrS_T:R
2. ttrR_T + R <--> ttrR_T:R        
        massaction: k_f(ttrR_T,R)=0.3*ttrR_T*R
        k_r(ttrR_T:R)=10*ttrR_T:R
3. ttrS + tt <--> ttrS_P + tt        
        massaction: k_f(ttrS,tt)=1.6*ttrS*tt
        k_r(ttrS_P,tt)=0.016*ttrS_P*tt
4. ttrS + ttrR <--> ttrR_ttrS        
        massaction: k_f(ttrS,ttrR)=0.0001*ttrS*ttrR
        k_r(ttrR_ttrS)=1*ttrR_ttrS
5. ttrS_P + ttrR <--> ttrR_ttrS_P        
        massaction: k_f(ttrS_P,ttrR)=0.0001*ttrS_P*ttrR
        k_r(ttrR_ttrS_P)=0.1*ttrR_ttrS_P
6. ttrR_ttrS_P <--> ttrR_P + ttrS        
        massaction: k_f(ttrR_ttrS_P)=1*ttrR_ttrS_P
        k_r(ttrR_P,ttrS)=1*ttrR_P*ttrS
7. 2 ttrR_P <--> ttrR:P_dimer        
        massaction: k_f(ttrR_P)=0.0083*ttrR_P^2
        k_r(ttrR:P_dimer)=5*ttrR:P_dimer
8. ttrR:P_dimer + pttr0 <--> pttr1        
        massaction: k_f(ttrR:P_dimer,pttr0)=0.3*ttrR:P_dimer*pttr0
        k_r(pttr1)=3*pttr1
9. pttr1_P --> GFP_T + P + pttr1        
        massaction: k_f(pttr1_P)=0.09*pttr1_P
10. GFP_T + R <--> GFP_T_R        
        massaction: k_f(GFP_T,R)=100*GFP_T*R
        k_r(GFP_T_R)=1*GFP_T_R
11. P:P1d --> ttrS_T + ttrR_T + P + P1d        
        massaction: k_f(P:P1d)=0.09*P:P1d
12. P + pttr1 <--> pttr1_P        
        massaction: k_f(P,pttr1)=0.1*P*pttr1
        k_r(pttr1_P)=1*pttr1_P
13. ttrS_T:R --> ttrS_T + ttrS + R        
        massaction: k_f(ttrS_T:R)=0.05*ttrS_T:R
14. ttrR_T:R --> ttrR_T + ttrR + R        
        massaction: k_f(ttrR_T:R)=0.05*ttrR_T:R
15. GFP_T_R --> GFP + GFP_T + R        
        massaction: k_f(GFP_T_R)=0.05*GFP_T_R
16. ttrS_T -->         
        massaction: k_f(ttrS_T)=0.9*ttrS_T
17. ttrR_T -->         
        massaction: k_f(ttrR_T)=0.9*ttrR_T
18. GFP_T -->         
        massaction: k_f(GFP_T)=0.9*GFP_T
19. ttrS -->         
        massaction: k_f(ttrS)=0.9*ttrS
20. ttrR -->         
        massaction: k_f(ttrR)=0.9*ttrR
21. GFP -->         
        massaction: k_f(GFP)=0.9*GFP
]
In [22]:
# An initial condition for each species (uninitialized species default to 0)
x0 = {
    "tt":200,
    "P1d":5,
    "pttr0":5,
    "P":25,
    "R":100,
}

timepoints = np.linspace(0, 600, 100)
In [23]:
Results = CRN_tt.simulate_with_bioscrape(timepoints,x0)
/Users/lianamerk/anaconda3/envs/python3.7-bioscrape/lib/python3.7/site-packages/biocrnpyler-0.1-py3.7.egg/biocrnpyler/chemical_reaction_network.py:870: UserWarning: The following species are uninitialized and their value has been defaulted to 0: ttrS, ttrR, 
In [24]:
Results.head()
Out[24]:
R P P1d P:P1d ttrS_T ttrR_T ttrR_T:R ttrS_T:R tt ttrS_P ... ttrR:P_dimer pttr0 pttr1 pttr1_P GFP_T GFP_T_R GFP ttrS ttrR time
0 100.000000 25.000000 5.000000 0.000000 0.000000 0.000000 0.000000 0.000000 200.0 0.000000 ... 0.000000e+00 5.0 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000 0.000000 0.000000
1 97.824944 20.002723 0.002723 4.997277 0.375266 0.375266 1.087528 1.087528 200.0 0.194098 ... 1.163234e-15 5.0 3.084548e-16 2.295029e-16 -1.763695e-18 -1.645280e-14 -9.145889e-16 0.002103 0.053622 6.060606
2 97.280408 20.002723 0.002723 4.997277 0.468971 0.468971 1.359796 1.359796 200.0 0.549631 ... 2.082295e-13 5.0 7.285205e-14 8.566282e-14 1.728165e-19 -5.795336e-15 -5.961676e-16 0.005691 0.073832 12.121212
3 97.146497 20.002723 0.002723 4.997277 0.492174 0.492174 1.426751 1.426751 200.0 0.928944 ... 2.437093e-12 5.0 9.621154e-13 1.344522e-12 4.207668e-17 2.741544e-13 9.917796e-15 0.009484 0.078838 18.181818
4 97.113649 20.002723 0.002723 4.997277 0.497875 0.497875 1.443175 1.443175 200.0 1.298786 ... 1.128285e-11 5.0 4.751947e-12 7.217950e-12 3.139122e-16 2.285017e-12 9.528891e-14 0.013174 0.080064 24.242424

5 rows × 23 columns

Holoviews plotting of Bioscrape results

Sort out Plotting

In [25]:
def plot_holoviews_ttr_bioscrape(results_melted):
    plot_group = []

    for row in results_melted.iterrows():
        sp = row[1]['species']

        if sp in ["P1d", "P:P1d", "pttr0", "pttr1", "pttr1_P"]:
            plot_group.append('promoters')

        elif sp in ["ttrS_T", "ttrR_T", "ttrR_T:R", "ttrS_T:R", "GFP_T", "GFP_T_R"]:
            plot_group.append('mRNA')

        elif sp in ["ttrS_P", "ttrR_ttrS", "ttrR_ttrS_P", "ttrR_P", "ttrR:P_dimer"]:
            plot_group.append('complexes')

        elif sp in ["ttrR", "ttrS", "GFP"]:
            plot_group.append('proteins')

        elif sp in ["P",   "R", "tt"]:
            plot_group.append('machinery and signals')
        else:
            print(sp)

    results_melted['plot_group'] = plot_group
    

    plot_list = []
    for group in results_melted['plot_group'].unique():
        plot_list.append(hv.Curve(
            data=results_melted.loc[results_melted['plot_group'] == group],
            kdims=['time', 'value'],
            vdims = ['species']
        ).groupby(
            [ 'species']
        ).opts(
            tools=['hover'],
#             legend_position = 'right',
            show_grid=True,
            title = f'{group}',
            width = 500,
            #shared_axes=False
        ).overlay('species').opts(legend_position = 'right',))
        
    return plot_list

Now, proceed with plotting

In [26]:
results_melted = pd.melt(Results, var_name='species', id_vars=['time'])
results_melted.head()
Out[26]:
time species value
0 0.000000 R 100.000000
1 6.060606 R 97.824944
2 12.121212 R 97.280408
3 18.181818 R 97.146497
4 24.242424 R 97.113649
In [27]:
plot_list = plot_holoviews_ttr_bioscrape(results_melted)
hv.Layout(plot_list).cols(2).opts(shared_axes=False)
Out[27]:

Tuning Function

In [28]:
def create_reactions(parameters):
    # P + P1d <-> P1d:P
    R1 = Reaction([P, P1d], [P_P1d], k = parameters['k1b'], k_rev = parameters['k1u'])

    # P1d:P -> ttrS_T + ttrR_T
    Rtx1 = Reaction([P_P1d], [ttrS_T, ttrR_T, P, P1d], k = parameters['ktx'])
    
    # ttrS_T + R <-> ttrS_T:R
    R2 = Reaction([ttrS_T, R], [ttrS_T_R], k = parameters['k2b'], k_rev = parameters['k2u'])

    # ttrS Translation
    Rtl1 = Reaction([ttrS_T_R], [ttrS_T, ttrS, R], k = parameters['ktl'])
    
    # ttrR_T + R <-> ttrR_T:R
    R3 = Reaction([ttrR_T, R], [ttrR_T_R], k = parameters['k3b'], k_rev = parameters['k3u'])

    # ttrR Translation
    Rtl2 = Reaction([ttrR_T_R], [ttrR_T, ttrR, R], k = parameters['ktl'])
    
    # ttrS + tt <-> ttrS_P + tt
    R4 = Reaction([ttrS, tt], [ttrS_P, tt], k = parameters['k4b'], k_rev = parameters['k4u'])
    
    # ttrR + ttrS <-> ttrR:ttrS
    R5 = Reaction([ttrS, ttrR], [ttrR_ttrS], k = parameters['k5b'], k_rev = parameters['k5u'])

    # ttrR + ttrS <-> ttrR:ttrS:P
    R6 = Reaction([ttrS_P, ttrR], [ttrR_ttrS_P], k = parameters['k6b'], k_rev = parameters['k6u'])

    # ttrR:ttrS:P <-> ttrR:P + ttrS
    R7 = Reaction([ttrR_ttrS_P], [ttrR_P, ttrS], k = parameters['k7b'], k_rev = parameters['k7u'])
    
    # ttrR:P -> ttrR
    rxn_dephos = Reaction([ttrR_P], [ttrR], k = parameters['k_dephos'])

    # 2ttrR:P <-> ttrR:P_dimer
    R8 = Reaction([ttrR_P, ttrR_P], [ttrR_P_dimer], k = parameters['k8b'], k_rev = parameters['k8u'])
    
    # ttrR:P_dimer + pttr0 <-> pttr1
    R9 = Reaction([ttrR_P_dimer, pttr0], [pttr1], k = parameters['k9b'], k_rev = parameters['k9u'])

    # P + pttr1 <-> pttr1:P
    Rtx2 = Reaction([P, pttr1], [pttr1_P], k = parameters['k10b'], k_rev = parameters['k10u'])

    # P + pttr1 <-> pttr1:P
    R10 = Reaction([pttr1_P], [GFP_T, P, pttr1], k = parameters['ktx'])

    # P + pttr1 <-> pttr1:P
    R11 = Reaction([GFP_T, R], [GFP_T_R], k = parameters['k11b'], k_rev = parameters['k11u'])

    # ttrS Translation
    Rtl3 = Reaction([GFP_T_R], [GFP, GFP_T, R], k = parameters['ktl'])

    # RNA deg
    rxnd_ttrS_T = Reaction([ttrS_T], [], k = parameters['delta'])
    rxnd_ttrR_T = Reaction([ttrR_T], [], k = parameters['delta'])
    rxnd_GFP_T = Reaction([GFP_T], [], k = parameters['delta'])

    # Protein Deg
    rxnd_ttrS = Reaction([ttrS], [], k = parameters['delta'])
    rxnd_ttrR = Reaction([ttrR], [], k = parameters['delta'])
    rxnd_GFP = Reaction([GFP], [], k = parameters['delta'])
    
    reactions = [R1, R4, R5, R6, R7, R8, R9, R10, R11,
                 Rtx1, Rtx2,
                 Rtl1, Rtl2, Rtl3,
                 rxnd_ttrS_T, rxnd_ttrR_T, rxnd_GFP_T,
                 rxnd_ttrS, rxnd_ttrR, rxnd_GFP]
    
    return reactions

RBS Tuning in Tetrathionate Circuit

We will change k2u, and k3u, the reaction rates that define ribosome unbinding from ttrS and ttrR transcripts, respectively. We will save the simulations as csv's to plot in another notebook such that the styles match well.

In [29]:
all_sims = pd.DataFrame()
In [30]:
#Parameters": List of tuples [("paramter name" (string)": value (float))]
parameters = {"k1b": 100, "k1u": 1, 
              "k2b": 3, "k2u": 10, 
              "k3b": 0.00003, "k3u": 10, 
              "k4b": 5, "k4u": 0.5,
              "k5b": 5, "k5u":6, 
              "k6b": 5, "k6u": 1, 
              "k7b": 0.01, "k7u": 1, 
              "k8b": 80, "k8u": 5, 
              "k9b": 0.004, "k9u": 3,
              "k10b": 50, "k10u":1, 
              "k11b": 10, "k11u": 1,     
              "ktx": 0.1, "ktl": 1,
              "k_dephos": 7, "delta": .9}
In [31]:
def plot_ttr_RBS_tuning(parameters, x0, k2u, k3u, title, timepoints):
    experimental_ttr_concs = np.array([ 0, 0.1, 0.25, 0.5, 0.75, 1, 5, 10])

    # Get color palette in order
    r = list(bokeh.palettes.viridis(10))
    r.reverse()
    r = r[2:]

    # ttrS_T + R <-> ttrS_T:R
    R2 = Reaction([ttrS_T, R], [ttrS_T_R], k = parameters['k2b'], k_rev = k2u)

    # ttrR_T + R <-> ttrR_T:R
    R3 = Reaction([ttrR_T, R], [ttrR_T_R], k = parameters['k3b'], k_rev = k3u)


    reactions = create_reactions(parameters)
    reactions.append(R2)
    reactions.append(R3)
    
    species = [R, P, P1d, P_P1d, ttrS_T, ttrR_T, ttrR_T_R, ttrS_T_R,
           tt, ttrS_P, ttrR_ttrS, ttrR_ttrS_P, ttrR_P, ttrR_P_dimer,
           pttr0, pttr1, pttr1_P, GFP_T, GFP_T_R, GFP]
    
    # Recompile the CRN with the right parameters
    CRN_curr = ChemicalReactionNetwork(species = species, reactions = reactions)
    
    p = bokeh.plotting.figure(width = 300, height= 250, x_axis_label='Time (hr)', y_axis_label = '[GFP]', title = title)
    i = 0

    for tt_conc in experimental_ttr_concs:
        x0["tt"] = tt_conc
        
        Results = CRN_curr.simulate_with_bioscrape(timepoints,x0)
        all_sims[f'{title}_{tt_conc}'] = Results['GFP']

        p.line(timepoints/60, Results['GFP'],color = r[i], legend_label = f'{tt_conc}', line_width = 2)
        i+=1

    p.legend.location = 'top_left'
    if k2u != 20:
        p.legend.visible = False
        
    p.y_range = Range1d(-0.1,2.5)

    return p
In [32]:
timepoints = np.linspace(0, 1400, 1400)
title_list = ['Very Low (LM15) Simulation', 'Low (LM14) Simulation', 'Medium (LM18) Simulation', 'High (LM19) Simulation']
In [33]:
x0 = {
    "tt":200,
    # PSC101 has copy number 5
    "P1d":5,
    "pttr0":5,
    # Polymerase: 15 uM assuming ~10000 Polymerase molecules / E. coli with a volume 1 um^3
    "P":15,
    # Ribosome: 120 uM assuming ~72000 Ribosomes / E. coli with a volume 1 um^3
    "R":120,
}
In [34]:
parameters = {"k1b": 10, "k1u": 0.0001, 
              "k2b": 0.3, "k2u": 10, 
              "k3b": 0.3, "k3u": 10, 
              "k4b": 1.6, "k4u": 0.016,
              "k5b": 0.0001, "k5u":1, 
              "k6b": 0.0001, "k6u": 0.1, 
              "k7b": 1, "k7u": 1, 
              "k8b": 0.0083, "k8u": 0.5, 
              "k9b": 0.3, "k9u": 3,
              "k10b": 10, "k10u":0.0001, 
              "k11b": 100, "k11u": 1,     
              "ktx": 0.1, "ktl":0.33,
              "k_dephos": 50, "delta": .9}
In [35]:
# Define higher off rates
k2u = [20, 12, 11, 10]
k3u = [15, 8, 6, 5]

plot_list = []

for i in range(len(k3u)):
    plot_list.append(plot_ttr_RBS_tuning(parameters, x0, k2u[i], k3u[i], title_list[i], timepoints))

bokeh.io.show(gridplot([plot_list]))
/Users/lianamerk/anaconda3/envs/python3.7-bioscrape/lib/python3.7/site-packages/biocrnpyler-0.1-py3.7.egg/biocrnpyler/chemical_reaction_network.py:870: UserWarning: The following species are uninitialized and their value has been defaulted to 0: ttrS, ttrR, 
In [36]:
all_sims['Time (hr)'] = timepoints
all_sims.to_csv('./data/all_sims.csv', index = False)

HRP AND Circuit

This $\sigma 54$-dependent split activator requires two proteins, hrpR and hrpS, to function. For now, let's just spike in the two proteins and see how our and gate responds.

In [37]:
txtl = CRNLab("GFP")

txtl.mixture("mixture1", extract = "TxTlExtract", mixture_volume = 1e-6, mixture_parameters = './parameters/BasicExtract.tsv')

#Define the activator and repressor species as proteins
prot1 = Protein("hrpR")
prot2 = Protein("hrpS")

dna = DNAassembly("comb_promoter",promoter=CombinatorialPromoter("phrpL",["hrpR","hrpS"], tx_capable_list = [['hrpS', 'hrpR'], ['hrpS'], ['hrpR']], leak=True),rbs="B0030",protein="GFP")

#TxTl Extract is a Mixture with more complex internal models
extract_1_TXTL = TxTlExtract(name = "e coli extract 1", components = [dna,prot1,prot2], parameter_file = "./parameters/parameters_and.txt")
CRN_extract_1 = extract_1_TXTL.compile_crn()

timepoints = np.linspace(0, 200, 1000)
x0 = {"dna_comb_promoter":5.0, "protein_hrpS":100, "protein_hrpR":100}
Re1 = CRN_extract_1.simulate_with_bioscrape(timepoints, initial_condition_dict = x0)

Re1.head()
Out[37]:
dna_comb_promoter complex_dna_comb_promoter_protein_hrpS complex_dna_comb_promoter_protein_RNAP complex_dna_comb_promoter_protein_hrpS_protein_RNAP complex_protein_RNAase_rna_comb_promoter complex_dna_comb_promoter_protein_hrpR_protein_hrpS_protein_RNAP protein_RNAase protein_GFP complex_dna_comb_promoter_protein_hrpR_protein_RNAP protein_Ribo complex_dna_comb_promoter_protein_hrpR_protein_hrpS complex_protein_Ribo_rna_comb_promoter complex_dna_comb_promoter_protein_hrpR rna_comb_promoter protein_RNAP protein_hrpR protein_hrpS time
0 5.000000 0.000000 0.000000 0.000000 0.000000 0.000000 6.000000 0.000000 0.000000 24.000000 0.000000 0.000000 0.000000 0.000000 3.000000 100.000000 100.000000 0.000000
1 0.000007 0.002111 0.009825 0.000002 0.004524 2.985622 5.995476 0.000466 0.000002 23.905593 2.000320 0.094407 0.002111 0.000595 0.004549 95.011945 95.011945 0.200200
2 0.000003 0.002106 0.001314 0.000002 0.013134 2.994570 5.986866 0.001879 0.000002 23.812736 1.999898 0.187264 0.002106 0.000981 0.004112 95.003424 95.003424 0.400400
3 0.000002 0.002105 0.000176 0.000002 0.025250 2.995767 5.974750 0.004204 0.000002 23.723280 1.999841 0.276720 0.002105 0.001357 0.004054 95.002285 95.002285 0.600601
4 0.000002 0.002105 0.000024 0.000002 0.040397 2.995927 5.959603 0.007409 0.000002 23.636834 1.999834 0.363166 0.002105 0.001724 0.004046 95.002133 95.002133 0.800801

Vary hrpR and hrpS proteins

Now, we will make sure the AND gate (combinatorial promoter) is working by varying protein concentrations.

Get plotting in order

Let's create a plotting function. We will plot a heatmap of protein GFP at the last slice of the simulation. Because we have no degradation of protein, the last timepoint will be the max. We will be able to pass in normalize, which divides all values by the maximum of all simulated values, or max_to_1, which will set the scale bar to be 0 to 1. This basically causes plots that might've been

In [38]:
def simulate_heatmap(inducer_1_name,
                     inducer_1_values,
                     inducer_2_name,
                     inducer_2_values,
                     title,
                     x0=x0,
                     normalize=True,
                     max_to_1 = False,
                     CRN=CRN_extract_1,
                     timepoints=timepoints):
    
    GFP_max = pd.DataFrame()
    
    for conc_1 in inducer_1_values:
        x0[inducer_1_name] = conc_1 
        
        for conc_2 in inducer_2_values:
            x0[inducer_2_name] = conc_2 #Change my initial condition dictionary
            
            Re1 = CRN_extract_1.simulate_with_bioscrape(timepoints, initial_condition_dict = x0)
            GFP_max = GFP_max.append({inducer_1_name:conc_1,
                            inducer_2_name:conc_2,
                            'GFP_max': Re1["protein_GFP"].values[-1]}, ignore_index=True)

    data = pd.pivot_table(data = GFP_max, index = inducer_1_name,
                   columns = inducer_2_name,
                   values = 'GFP_max')
    
    fig, ax = plt.subplots(figsize=(6, 4.5))

    if normalize:
        data /= np.max(GFP_max['GFP_max'])
        
    if max_to_1:
        heat_map = sns.heatmap(data,
                         cmap = 'viridis',
                        linewidths=.1, vmin = 0, vmax = 1,
                        cbar_kws={'label': 'Normalized GFP'})
    else:
        heat_map = sns.heatmap(data,
                         cmap = 'viridis',
                        linewidths=.2,
                        cbar_kws={'label': 'Normalized GFP'})
        
        
    sns.set(font_scale=1.4)
    # plt.tick_params(axis='both', labelsize=20)

    plt.xticks(fontsize=18, fontstyle='normal', rotation = 90)
    plt.yticks(fontsize=18, fontstyle='normal')


    plt.xlabel(inducer_1_name, fontsize=18, fontweight='bold')
    plt.ylabel(inducer_2_name, fontsize=18, fontweight='bold')


    heat_map.set_title(title)
    fig.tight_layout()
    
    return fig, heat_map

And Gate with leak in hrpR

Let's say hrpR mutated and now it can recruit RNAP better than it did beter. We will change phrpL_hrpR_RNAP ku = 5 instead of 500.

In [39]:
hrpS_values = [0.0001, 0.001, 0.01, 0.1, 0.5, 1]
hrpR_values = [0.0001, 0.001, 0.01, 0.1, 0.5, 1]

dna = DNAassembly("comb_promoter",promoter=CombinatorialPromoter("phrpL",["hrpR","hrpS"], tx_capable_list = [['hrpS', 'hrpR'], ['hrpS'], ['hrpR']], leak=True),rbs="B0030",protein="GFP")
extract_1_TXTL = TxTlExtract(name = "e coli extract 1", components = [dna,prot1,prot2], parameter_file = "./parameters/parameters_and_leak5_hrpR.txt")
CRN_extract_1 = extract_1_TXTL.compile_crn()

timepoints = np.linspace(0, 200, 1000)
x0 = {"dna_comb_promoter":5.0, "protein_hrpS":100, "protein_hrpR":100}

fig, heat_map = simulate_heatmap('protein_hrpS',hrpS_values,'protein_hrpR',hrpR_values,x0=x0,normalize=True,title='Simulated Protein-Based \n AND gate 1x $k_u^{hrpR}$', CRN=CRN_extract_1)
heat_map.figure.axes[-1].yaxis.label.set_size(20)

heat_map.set(xlabel='hrpS protein', ylabel='hrpR protein')
heat_map.invert_yaxis()
fig.savefig("./figures/fig5/fig5c_1.png", bbox_inches = "tight")
plt.show()
In [40]:
hrpS_values = [0.0001, 0.001, 0.01, 0.1, 0.5, 1]
hrpR_values = [0.0001, 0.001, 0.01, 0.1, 0.5, 1]

dna = DNAassembly("comb_promoter",promoter=CombinatorialPromoter("phrpL",["hrpR","hrpS"], tx_capable_list = [['hrpS', 'hrpR'], ['hrpS'], ['hrpR']], leak=True),rbs="B0030",protein="GFP")
extract_1_TXTL = TxTlExtract(name = "e coli extract 1", components = [dna,prot1,prot2], parameter_file = "./parameters/parameters_and_leak50_hrpR.txt")
CRN_extract_1 = extract_1_TXTL.compile_crn()

timepoints = np.linspace(0, 200, 1000)
x0 = {"dna_comb_promoter":5.0, "protein_hrpS":100, "protein_hrpR":100}

fig, heat_map = simulate_heatmap('protein_hrpS',hrpS_values,'protein_hrpR',hrpR_values,x0=x0,normalize=True,title='Simulated Protein-Based \n AND gate 10x $k_u^{hrpR}$', CRN=CRN_extract_1)
heat_map.figure.axes[-1].yaxis.label.set_size(20)

heat_map.set(xlabel='hrpS protein', ylabel='hrpR protein')
heat_map.invert_yaxis()
fig.savefig("./figures/fig5/fig5c_2.png", bbox_inches = "tight")
plt.show()
In [41]:
hrpS_values = [0.0001, 0.001, 0.01, 0.1, 0.5, 1]
hrpR_values = [0.0001, 0.001, 0.01, 0.1, 0.5, 1]

dna = DNAassembly("comb_promoter",promoter=CombinatorialPromoter("phrpL",["hrpR","hrpS"], tx_capable_list = [['hrpS', 'hrpR'], ['hrpS'], ['hrpR']], leak=True),rbs="B0030",protein="GFP")
extract_1_TXTL = TxTlExtract(name = "e coli extract 1", components = [dna,prot1,prot2], parameter_file = "./parameters/parameters_and.txt")
CRN_extract_1 = extract_1_TXTL.compile_crn()

timepoints = np.linspace(0, 200, 1000)
x0 = {"dna_comb_promoter":5.0, "protein_hrpS":100, "protein_hrpR":100}

fig, heat_map = simulate_heatmap('protein_hrpS',hrpS_values,'protein_hrpR',hrpR_values,x0=x0,normalize=True,title='Simulated Protein-Based \n AND gate 100x $k_u^{hrpR}$', CRN=CRN_extract_1)
heat_map.figure.axes[-1].yaxis.label.set_size(20)

heat_map.set(xlabel='hrpS protein', ylabel='hrpR protein')
heat_map.invert_yaxis()
fig.savefig("./figures/fig5/fig5c_3.png", bbox_inches = "tight")
plt.show()

Full Compiled AND gate with expression of hrpR, S

LacI promoter class

In [42]:
class LacIPromoter2(Promoter):
    def __init__(self, name = "pLac",  IPTG = "IPTG",
                 lacI = "lacI", leak = True, assembly = None,
                 transcript = None, length = 0, mechanisms = {},
                 parameters = {}, **keywords):

        Promoter.__init__(self, name=name, transcript=transcript, assembly=assembly, length=length, parameters=parameters, **keywords)
        #RegulatedPromoter.__init__(self = self, name = name, regulators = regulators)

        # Transcription is already included in promoter class
        self.default_mechanisms = {"binding": One_Step_Cooperative_Binding()}

        self.leak = leak

        self.IPTG = self.set_species(IPTG, material_type = "inducer")
        self.lacI = self.set_species(lacI, material_type = "protein")

        self.plac_2_lacI = None
        self.IPTG_2_lacI = None

    def update_species(self):
        mech_tx = self.mechanisms['transcription']
        mech_b = self.mechanisms['binding']
        species = []
        self.complexes = []

        """A reaction where n binders (s1) bind to 1 bindee (s2) in two steps
           2lacI <--> 2xlacI
           2xlacI + pLac <--> 2xlacI:pLac
        """
        species_b = mech_b.update_species(self.lacI, self.assembly.dna, n_mer_species=self.lacI, component = self, part_id = self.assembly.dna.name+'_2x_'+self.lacI.name, cooperativity=2)
        species += species_b

        #print('plac stuff \n',species_b)
        for s in species_b:
            if isinstance(s, ComplexSpecies) and s.species.count(self.lacI) == 2 and self.assembly.dna in s.species:
                self.plac_2_lacI = s

        if self.leak is not False:
            species += mech_tx.update_species(dna = self.plac_2_lacI, component = self, part_id = self.name, \
                                        transcript = self.transcript, protein = self.assembly.protein)

        # IPTG can bind to lacI to sequester it
        # give citation for 1 being enough
        """A reaction where n binders (s1) bind to 1 bindee (s2) in two steps
           2lacI <--> 2xlacI
           2xlacI + IPTG <--> 2xlacI:IPTG
        """
        species_b = mech_b.update_species(self.lacI, self.IPTG, component = self, part_id = self.IPTG.name+'_2x_'+self.lacI.name, cooperativity=2)
        #print('iptg stuff \n', species_b)

        species += species_b

        for s in species_b:
            if isinstance(s, ComplexSpecies) and s.species.count(self.lacI) == 2 and self.IPTG in s.species:
                self.IPTG_2_lacI = s

        # Unbound transcription
        species += mech_tx.update_species(dna = self.assembly.dna, component = self, part_id = self.name, \
                                        transcript = self.transcript, protein = self.assembly.protein)

        return species


    def update_reactions(self):
        print('r')
        mech_tx = self.mechanisms['transcription']
        mech_b = self.mechanisms['binding']

        if self.IPTG_2_lacI == None or self.plac_2_lacI == None:
            self.update_species()

        if self.leak != False:
            # Will need transcription, plac_2x_lacI, ktx, ku, kb
            reactions += mech_tx.update_reactions(dna = self.plac_2_lacI, component = self, part_id = self.plac_2_lacI.name, \
                                        transcript = self.transcript, protein = self.assembly.protein)

        reactions = []

        # Repression
        # Will need two_step_cooperative_binding, plac_2x_lacI, kb1 ku1, kb2, ku2
        reactions += mech_b.update_reactions(self.lacI, self.assembly.dna, component = self, \
                                                                part_id = self.plac_2_lacI.name, cooperativity = 2)


        # LacI binding to IPTG
        # Will need two_step_cooperative_binding, IPTG_2x_lacI, kb1 ku1, kb2, ku2
        # .name gets rid of the complex
        # .type will give you complex, species, etc
        # .attributes gives you properties like phosphorylation state
        reactions += mech_b.update_reactions(self.lacI, self.IPTG, component = self, part_id = self.IPTG_2_lacI.name, cooperativity=2)

        '''
        Any bit of IPTG destroys the bound structure
        IPTG + plac_2xlacI --> IPTG:2xlacI + plac
        '''
        if isinstance(mech_b, Two_Step_Cooperative_Binding):
            kb = self.get_parameter("kb2", part_id = self.IPTG_2_lacI.name, mechanism = mech_b)
        else:
            kb = self.get_parameter("kb", part_id = self.IPTG_2_lacI.name, mechanism = mech_b)
            print(self.IPTG_2_lacI.name, kb)

        reactions.append(Reaction([self.IPTG, self.plac_2_lacI], [self.IPTG_2_lacI, self.assembly.dna], k=kb))

        # Unbound transcription plac, kb, ku, ktx
        reactions += mech_tx.update_reactions(dna = self.assembly.dna, component = self, part_id = self.name, \
                                        transcript = self.transcript, protein = self.assembly.protein)

        return reactions

Sort out parameter printer

In [43]:
def param_printer2(x0, parameters):
    inits = pd.DataFrame(list(x0.items()),columns = ['param','value']) 
    

    df = pd.DataFrame(list(parameters.items()),columns = ['param','value']) 
    col_part_id = []
    mech = []
    for val in df['param']:
        if isinstance(val, tuple):
            col_part_id.append(val[0])
            mech.append(val[1])
        else:
            col_part_id.append(' ')
            mech.append(val)

    df['param'] = mech
    df['part'] = col_part_id
    df.reindex(['part','param','value'],axis=1)
    df[''] = ['|']*len(col_part_id)
    
    return pd.concat([df, inits],  axis=1).fillna("")

AND simulation

In [44]:
parameters = {
    ('pttr_phosphate_P_protein_ttrR', 'cooperativity'): 2.0, 
    ('pttr_phosphate_P_protein_ttrR', 'ktx'): 0.17025,  #this is the transcription rate
    ('pttr_phosphate_P_protein_ttrR', 'ku'): 11.01321586, #this is the polymerase unbinding rate

    ('transcription_mm',  'ktx'): .01,    #These are the parameters for transcription leak
    ('transcription_mm','kb'):100, #default binding rate!
    #('ttrS-P',"kdephos"):1000, #auto dephosphorylation
    ('ligand_tt_protein_ttrS-P','kdephos'):1000, #dephosphorylation rate of ttrS
    #('ttrS-P',"kphos"):10, #ttrs phosphorylates itself automatically. This is a proxy for sensing
    ('ligand_tt_protein_ttrS-P','kphos'):10, #equivalent to kphos, from before
    #('ttrS-ttrR','ktransfer'):5, #rate of transfer from ttrS to ttrR
    ('ligand_tt_protein_ttrS-ttrR', 'ktransfer'):5, #this is the new transfer rate
    ('ttrR-P',"kdephos"):1000, #auto dephosphorylation
    "kb":100, "ku":10, "ktx":.05, "ktl":.2, "kdeg":3, 'cooperativity':2, #some default parameters
    
    ('phrpL_hrpR', 'cooperativity'): 1.0, #default is 2, but these proteins bind at cooperativity = 1
    ('phrpL_hrpS', 'cooperativity'): 1.0, 
    ('phrpL_hrpR_hrpS_RNAP', 'ktx'): 0.17025,  #this is the transcription rate
    ('phrpL_hrpR_hrpS_RNAP', 'ku'): 11.01321586, #this is the polymerase unbinding rate

    ('translation_mm', 'B0030', 'ku'): 888, #Unbinding
    ('translation_mm', 'B0030', 'ktl'): .15, #Translation Rate

    ('transcription_mm',  'ktx'): .01,    #These are the parameters for transcription leak
    ('transcription_mm', 'ku'): 100,
    
    ('phrpL_hrpR','ku'):100,
    ('phrpL_hrpS','ku'):100,
    ('phrpL_hrpS_RNAP', 'ku'): 1000.,#we'll say that leak has 100x the unbinding rate
    ('phrpL_hrpR_RNAP', 'ku'): 1000,
    ('phrpL','ku'):100000,
    
    ('binding', 'pLac_lacI','ktx'):0.001,
    ('plac_lacI','kb'):10,
    ('plac_lacI','ku'):0.001,
    
    ('binding', 'plac_hrpS_2x_lacI','kb'):8888,
    ("simple_transcription", 'plac','ktx'):0.1,
    
    #Make a lot of repressor
    ('translation_mm', 'rbsL', 'ku'): 0.0001, #Unbinding
    ('translation_mm', 'rbsL', 'ktl'): 2, #Translation Rate
    ('plac_hrpS','kb'):100,
    ('plac_hrpS','ku'):0.0001
}


#P_rbsS_ttrS
constit_ttrS = DNAassembly(name = "constit_ttrS", promoter = "P", rbs = "rbsS", protein = "ttrS")
#P_rbsR_ttrR
constit_ttrR = DNAassembly(name = "constit_ttrR", promoter = "P", rbs = "rbsR", protein = "ttrR")
#P_arl_lacI
constit_lacI = DNAassembly(name = "constit_lacI", promoter = "P", rbs = "rbsL", protein = "lacI")
#P_arl_lacI
constit_mScarlet = DNAassembly(name = "constit_mScarlet", promoter = "P", rbs = "B0034", protein = "mScarlet")

# plac_hrpS
IPTG = Species("IPTG", material_type='inducer')
#IPTG = Protein("IPTG")

lac_p = LacIPromoter2(name="pLac", IPTG = "IPTG", lacI = "lacI", leak=False)
plac_hrpS = DNAassembly(name = "plac_hrpS", promoter=lac_p, rbs="rbsH",protein="hrpS")

# pttr_hrpR
tt = Species("tt",material_type="ligand")
ttrSbound = ChemicalComplex([Protein("ttrS"),tt])
ttrS3 = PhosphoTransferase(ttrSbound.species,targets=["ttrR"])
ttrR = PhosphoProtein("ttrR")
# Need to change from UTR1
pttr_hrpR = DNAassembly(name = "pttr_hrpR",promoter=RegulatedPromoter("pttr",[ttrR.get_phosphorylated_species()],leak=False),rbs="UTR1",protein="hrpR")

hrpL_gfp = DNAassembly("hrpL_gfp",promoter=CombinatorialPromoter("phrpL",["hrpR","hrpS"], tx_capable_list = [['hrpS', 'hrpR'], ['hrpS'], ['hrpR']], leak=False),rbs="B0030",protein="GFP")

extract_1_TXTL = TxTlExtract(name = "e coli", components = [pttr_hrpR,ttrSbound,ttrS3,ttrR,constit_ttrS,constit_ttrR,
                                                                     constit_lacI, plac_hrpS, hrpL_gfp, constit_mScarlet],parameters=parameters)

CRN_extract_1 = extract_1_TXTL.compile_crn()
r
inducer_IPTG_2x_protein_lacI 100
In [45]:
param_printer2(x0, parameters)
Out[45]:
param value part param value
0 cooperativity 2.000000 pttr_phosphate_P_protein_ttrR | dna_comb_promoter 5
1 ktx 0.170250 pttr_phosphate_P_protein_ttrR | protein_hrpS 1
2 ku 11.013216 pttr_phosphate_P_protein_ttrR | protein_hrpR 1
3 ktx 0.010000 transcription_mm |
4 kb 100.000000 transcription_mm |
5 kdephos 1000.000000 ligand_tt_protein_ttrS-P |
6 kphos 10.000000 ligand_tt_protein_ttrS-P |
7 ktransfer 5.000000 ligand_tt_protein_ttrS-ttrR |
8 kdephos 1000.000000 ttrR-P |
9 kb 100.000000 |
10 ku 10.000000 |
11 ktx 0.050000 |
12 ktl 0.200000 |
13 kdeg 3.000000 |
14 cooperativity 2.000000 |
15 cooperativity 1.000000 phrpL_hrpR |
16 cooperativity 1.000000 phrpL_hrpS |
17 ktx 0.170250 phrpL_hrpR_hrpS_RNAP |
18 ku 11.013216 phrpL_hrpR_hrpS_RNAP |
19 B0030 888.000000 translation_mm |
20 B0030 0.150000 translation_mm |
21 ku 100.000000 transcription_mm |
22 ku 100.000000 phrpL_hrpR |
23 ku 100.000000 phrpL_hrpS |
24 ku 1000.000000 phrpL_hrpS_RNAP |
25 ku 1000.000000 phrpL_hrpR_RNAP |
26 ku 100000.000000 phrpL |
27 pLac_lacI 0.001000 binding |
28 kb 10.000000 plac_lacI |
29 ku 0.001000 plac_lacI |
30 plac_hrpS_2x_lacI 8888.000000 binding |
31 plac 0.100000 simple_transcription |
32 rbsL 0.000100 translation_mm |
33 rbsL 2.000000 translation_mm |
34 kb 100.000000 plac_hrpS |
35 ku 0.000100 plac_hrpS |
In [46]:
x0 = {"dna_pttr_hrpR":1,
      "dna_constit_ttrS":1,
      "dna_constit_ttrR":1,
      "dna_constit_lacI":1,
      "dna_plac_hrpS":1,
      "dna_hrpL_gfp":1,
      "dna_constit_mScarlet":5,
      "phosphate_P":100,
      "protein_RNAP":15,
      "protein_RNAase":30,
      "protein_Ribo":120,
     "protein_lacI": 5}
In [47]:
IPTG_um_list = np.logspace(-3, 2, 6).round(decimals=6)
tt_um_list = np.logspace(-3, 2, 6).round(decimals=6)

fig, heat_map = simulate_heatmap('ligand_tt',tt_um_list,'inducer_IPTG',IPTG_um_list,x0=x0,normalize=True,title='Compiled Full AND gate \n lacI Spike', CRN=CRN_extract_1)

heat_map.set(xlabel='IPTG (uM)', ylabel='Tetrathionate (uM)')
heat_map.invert_yaxis()

plt.yticks(fontsize=18, fontstyle='normal', rotation = 0)
fig.savefig("./figures/fig6/fig6d.png")
plt.show()
/Users/lianamerk/anaconda3/envs/python3.7-bioscrape/lib/python3.7/site-packages/biocrnpyler-0.1-py3.7.egg/biocrnpyler/chemical_reaction_network.py:870: UserWarning: The following species are uninitialized and their value has been defaulted to 0: inducer_IPTG, 
/Users/lianamerk/anaconda3/envs/python3.7-bioscrape/lib/python3.7/site-packages/scipy/integrate/odepack.py:248: ODEintWarning: Excess work done on this call (perhaps wrong Dfun type). Run with full_output = 1 to get quantitative information.
  warnings.warn(warning_msg, ODEintWarning)
odeint failed with mxstep=500...odeint failed with mxstep=500...odeint failed with mxstep=500...odeint failed with mxstep=500...odeint failed with mxstep=500...odeint failed with mxstep=500...
100) { clearInterval(timer); console.log("Bokeh: ERROR: Unable to run BokehJS code because BokehJS library is missing"); } } }, 10, root) } })(window);
In [36]:
all_sims['Time (hr)'] = timepoints
all_sims.to_csv('./data/all_sims.csv', index = False)

HRP AND Circuit

This $\sigma 54$-dependent split activator requires two proteins, hrpR and hrpS, to function. For now, let's just spike in the two proteins and see how our and gate responds.

In [37]:
txtl = CRNLab("GFP")

txtl.mixture("mixture1", extract = "TxTlExtract", mixture_volume = 1e-6, mixture_parameters = './parameters/BasicExtract.tsv')

#Define the activator and repressor species as proteins
prot1 = Protein("hrpR")
prot2 = Protein("hrpS")

dna = DNAassembly("comb_promoter",promoter=CombinatorialPromoter("phrpL",["hrpR","hrpS"], tx_capable_list = [['hrpS', 'hrpR'], ['hrpS'], ['hrpR']], leak=True),rbs="B0030",protein="GFP")

#TxTl Extract is a Mixture with more complex internal models
extract_1_TXTL = TxTlExtract(name = "e coli extract 1", components = [dna,prot1,prot2], parameter_file = "./parameters/parameters_and.txt")
CRN_extract_1 = extract_1_TXTL.compile_crn()

timepoints = np.linspace(0, 200, 1000)
x0 = {"dna_comb_promoter":5.0, "protein_hrpS":100, "protein_hrpR":100}
Re1 = CRN_extract_1.simulate_with_bioscrape(timepoints, initial_condition_dict = x0)

Re1.head()
Out[37]:
dna_comb_promoter complex_dna_comb_promoter_protein_hrpS complex_dna_comb_promoter_protein_RNAP complex_dna_comb_promoter_protein_hrpS_protein_RNAP complex_protein_RNAase_rna_comb_promoter complex_dna_comb_promoter_protein_hrpR_protein_hrpS_protein_RNAP protein_RNAase protein_GFP complex_dna_comb_promoter_protein_hrpR_protein_RNAP protein_Ribo complex_dna_comb_promoter_protein_hrpR_protein_hrpS complex_protein_Ribo_rna_comb_promoter complex_dna_comb_promoter_protein_hrpR rna_comb_promoter protein_RNAP protein_hrpR protein_hrpS time
0 5.000000 0.000000 0.000000 0.000000 0.000000 0.000000 6.000000 0.000000 0.000000 24.000000 0.000000 0.000000 0.000000 0.000000 3.000000 100.000000 100.000000 0.000000
1 0.000007 0.002111 0.009825 0.000002 0.004524 2.985622 5.995476 0.000466 0.000002 23.905593 2.000320 0.094407 0.002111 0.000595 0.004549 95.011945 95.011945 0.200200
2 0.000003 0.002106 0.001314 0.000002 0.013134 2.994570 5.986866 0.001879 0.000002 23.812736 1.999898 0.187264 0.002106 0.000981 0.004112 95.003424 95.003424 0.400400
3 0.000002 0.002105 0.000176 0.000002 0.025250 2.995767 5.974750 0.004204 0.000002 23.723280 1.999841 0.276720 0.002105 0.001357 0.004054 95.002285 95.002285 0.600601
4 0.000002 0.002105 0.000024 0.000002 0.040397 2.995927 5.959603 0.007409 0.000002 23.636834 1.999834 0.363166 0.002105 0.001724 0.004046 95.002133 95.002133 0.800801

Vary hrpR and hrpS proteins

Now, we will make sure the AND gate (combinatorial promoter) is working by varying protein concentrations.

Get plotting in order

Let's create a plotting function. We will plot a heatmap of protein GFP at the last slice of the simulation. Because we have no degradation of protein, the last timepoint will be the max. We will be able to pass in normalize, which divides all values by the maximum of all simulated values, or max_to_1, which will set the scale bar to be 0 to 1. This basically causes plots that might've been

In [38]:
def simulate_heatmap(inducer_1_name,
                     inducer_1_values,
                     inducer_2_name,
                     inducer_2_values,
                     title,
                     x0=x0,
                     normalize=True,
                     max_to_1 = False,
                     CRN=CRN_extract_1,
                     timepoints=timepoints):
    
    GFP_max = pd.DataFrame()
    
    for conc_1 in inducer_1_values:
        x0[inducer_1_name] = conc_1 
        
        for conc_2 in inducer_2_values:
            x0[inducer_2_name] = conc_2 #Change my initial condition dictionary
            
            Re1 = CRN_extract_1.simulate_with_bioscrape(timepoints, initial_condition_dict = x0)
            GFP_max = GFP_max.append({inducer_1_name:conc_1,
                            inducer_2_name:conc_2,
                            'GFP_max': Re1["protein_GFP"].values[-1]}, ignore_index=True)

    data = pd.pivot_table(data = GFP_max, index = inducer_1_name,
                   columns = inducer_2_name,
                   values = 'GFP_max')
    
    fig, ax = plt.subplots(figsize=(6, 4.5))

    if normalize:
        data /= np.max(GFP_max['GFP_max'])
        
    if max_to_1:
        heat_map = sns.heatmap(data,
                         cmap = 'viridis',
                        linewidths=.1, vmin = 0, vmax = 1,
                        cbar_kws={'label': 'Normalized GFP'})
    else:
        heat_map = sns.heatmap(data,
                         cmap = 'viridis',
                        linewidths=.2,
                        cbar_kws={'label': 'Normalized GFP'})
        
        
    sns.set(font_scale=1.4)
    # plt.tick_params(axis='both', labelsize=20)

    plt.xticks(fontsize=18, fontstyle='normal', rotation = 90)
    plt.yticks(fontsize=18, fontstyle='normal')


    plt.xlabel(inducer_1_name, fontsize=18, fontweight='bold')
    plt.ylabel(inducer_2_name, fontsize=18, fontweight='bold')


    heat_map.set_title(title)
    fig.tight_layout()
    
    return fig, heat_map

And Gate with leak in hrpR

Let's say hrpR mutated and now it can recruit RNAP better than it did beter. We will change phrpL_hrpR_RNAP ku = 5 instead of 500.

In [39]:
hrpS_values = [0.0001, 0.001, 0.01, 0.1, 0.5, 1]
hrpR_values = [0.0001, 0.001, 0.01, 0.1, 0.5, 1]

dna = DNAassembly("comb_promoter",promoter=CombinatorialPromoter("phrpL",["hrpR","hrpS"], tx_capable_list = [['hrpS', 'hrpR'], ['hrpS'], ['hrpR']], leak=True),rbs="B0030",protein="GFP")
extract_1_TXTL = TxTlExtract(name = "e coli extract 1", components = [dna,prot1,prot2], parameter_file = "./parameters/parameters_and_leak5_hrpR.txt")
CRN_extract_1 = extract_1_TXTL.compile_crn()

timepoints = np.linspace(0, 200, 1000)
x0 = {"dna_comb_promoter":5.0, "protein_hrpS":100, "protein_hrpR":100}

fig, heat_map = simulate_heatmap('protein_hrpS',hrpS_values,'protein_hrpR',hrpR_values,x0=x0,normalize=True,title='Simulated Protein-Based \n AND gate 1x $k_u^{hrpR}$', CRN=CRN_extract_1)
heat_map.figure.axes[-1].yaxis.label.set_size(20)

heat_map.set(xlabel='hrpS protein', ylabel='hrpR protein')
heat_map.invert_yaxis()
fig.savefig("./figures/fig5/fig5c_1.png", bbox_inches = "tight")
plt.show()
In [40]:
hrpS_values = [0.0001, 0.001, 0.01, 0.1, 0.5, 1]
hrpR_values = [0.0001, 0.001, 0.01, 0.1, 0.5, 1]

dna = DNAassembly("comb_promoter",promoter=CombinatorialPromoter("phrpL",["hrpR","hrpS"], tx_capable_list = [['hrpS', 'hrpR'], ['hrpS'], ['hrpR']], leak=True),rbs="B0030",protein="GFP")
extract_1_TXTL = TxTlExtract(name = "e coli extract 1", components = [dna,prot1,prot2], parameter_file = "./parameters/parameters_and_leak50_hrpR.txt")
CRN_extract_1 = extract_1_TXTL.compile_crn()

timepoints = np.linspace(0, 200, 1000)
x0 = {"dna_comb_promoter":5.0, "protein_hrpS":100, "protein_hrpR":100}

fig, heat_map = simulate_heatmap('protein_hrpS',hrpS_values,'protein_hrpR',hrpR_values,x0=x0,normalize=True,title='Simulated Protein-Based \n AND gate 10x $k_u^{hrpR}$', CRN=CRN_extract_1)
heat_map.figure.axes[-1].yaxis.label.set_size(20)

heat_map.set(xlabel='hrpS protein', ylabel='hrpR protein')
heat_map.invert_yaxis()
fig.savefig("./figures/fig5/fig5c_2.png", bbox_inches = "tight")
plt.show()
In [41]:
hrpS_values = [0.0001, 0.001, 0.01, 0.1, 0.5, 1]
hrpR_values = [0.0001, 0.001, 0.01, 0.1, 0.5, 1]

dna = DNAassembly("comb_promoter",promoter=CombinatorialPromoter("phrpL",["hrpR","hrpS"], tx_capable_list = [['hrpS', 'hrpR'], ['hrpS'], ['hrpR']], leak=True),rbs="B0030",protein="GFP")
extract_1_TXTL = TxTlExtract(name = "e coli extract 1", components = [dna,prot1,prot2], parameter_file = "./parameters/parameters_and.txt")
CRN_extract_1 = extract_1_TXTL.compile_crn()

timepoints = np.linspace(0, 200, 1000)
x0 = {"dna_comb_promoter":5.0, "protein_hrpS":100, "protein_hrpR":100}

fig, heat_map = simulate_heatmap('protein_hrpS',hrpS_values,'protein_hrpR',hrpR_values,x0=x0,normalize=True,title='Simulated Protein-Based \n AND gate 100x $k_u^{hrpR}$', CRN=CRN_extract_1)
heat_map.figure.axes[-1].yaxis.label.set_size(20)

heat_map.set(xlabel='hrpS protein', ylabel='hrpR protein')
heat_map.invert_yaxis()
fig.savefig("./figures/fig5/fig5c_3.png", bbox_inches = "tight")
plt.show()

Full Compiled AND gate with expression of hrpR, S

LacI promoter class

In [42]:
class LacIPromoter2(Promoter):
    def __init__(self, name = "pLac",  IPTG = "IPTG",
                 lacI = "lacI", leak = True, assembly = None,
                 transcript = None, length = 0, mechanisms = {},
                 parameters = {}, **keywords):

        Promoter.__init__(self, name=name, transcript=transcript, assembly=assembly, length=length, parameters=parameters, **keywords)
        #RegulatedPromoter.__init__(self = self, name = name, regulators = regulators)

        # Transcription is already included in promoter class
        self.default_mechanisms = {"binding": One_Step_Cooperative_Binding()}

        self.leak = leak

        self.IPTG = self.set_species(IPTG, material_type = "inducer")
        self.lacI = self.set_species(lacI, material_type = "protein")

        self.plac_2_lacI = None
        self.IPTG_2_lacI = None

    def update_species(self):
        mech_tx = self.mechanisms['transcription']
        mech_b = self.mechanisms['binding']
        species = []
        self.complexes = []

        """A reaction where n binders (s1) bind to 1 bindee (s2) in two steps
           2lacI <--> 2xlacI
           2xlacI + pLac <--> 2xlacI:pLac
        """
        species_b = mech_b.update_species(self.lacI, self.assembly.dna, n_mer_species=self.lacI, component = self, part_id = self.assembly.dna.name+'_2x_'+self.lacI.name, cooperativity=2)
        species += species_b

        #print('plac stuff \n',species_b)
        for s in species_b:
            if isinstance(s, ComplexSpecies) and s.species.count(self.lacI) == 2 and self.assembly.dna in s.species:
                self.plac_2_lacI = s

        if self.leak is not False:
            species += mech_tx.update_species(dna = self.plac_2_lacI, component = self, part_id = self.name, \
                                        transcript = self.transcript, protein = self.assembly.protein)

        # IPTG can bind to lacI to sequester it
        # give citation for 1 being enough
        """A reaction where n binders (s1) bind to 1 bindee (s2) in two steps
           2lacI <--> 2xlacI
           2xlacI + IPTG <--> 2xlacI:IPTG
        """
        species_b = mech_b.update_species(self.lacI, self.IPTG, component = self, part_id = self.IPTG.name+'_2x_'+self.lacI.name, cooperativity=2)
        #print('iptg stuff \n', species_b)

        species += species_b

        for s in species_b:
            if isinstance(s, ComplexSpecies) and s.species.count(self.lacI) == 2 and self.IPTG in s.species:
                self.IPTG_2_lacI = s

        # Unbound transcription
        species += mech_tx.update_species(dna = self.assembly.dna, component = self, part_id = self.name, \
                                        transcript = self.transcript, protein = self.assembly.protein)

        return species


    def update_reactions(self):
        print('r')
        mech_tx = self.mechanisms['transcription']
        mech_b = self.mechanisms['binding']

        if self.IPTG_2_lacI == None or self.plac_2_lacI == None:
            self.update_species()

        if self.leak != False:
            # Will need transcription, plac_2x_lacI, ktx, ku, kb
            reactions += mech_tx.update_reactions(dna = self.plac_2_lacI, component = self, part_id = self.plac_2_lacI.name, \
                                        transcript = self.transcript, protein = self.assembly.protein)

        reactions = []

        # Repression
        # Will need two_step_cooperative_binding, plac_2x_lacI, kb1 ku1, kb2, ku2
        reactions += mech_b.update_reactions(self.lacI, self.assembly.dna, component = self, \
                                                                part_id = self.plac_2_lacI.name, cooperativity = 2)


        # LacI binding to IPTG
        # Will need two_step_cooperative_binding, IPTG_2x_lacI, kb1 ku1, kb2, ku2
        # .name gets rid of the complex
        # .type will give you complex, species, etc
        # .attributes gives you properties like phosphorylation state
        reactions += mech_b.update_reactions(self.lacI, self.IPTG, component = self, part_id = self.IPTG_2_lacI.name, cooperativity=2)

        '''
        Any bit of IPTG destroys the bound structure
        IPTG + plac_2xlacI --> IPTG:2xlacI + plac
        '''
        if isinstance(mech_b, Two_Step_Cooperative_Binding):
            kb = self.get_parameter("kb2", part_id = self.IPTG_2_lacI.name, mechanism = mech_b)
        else:
            kb = self.get_parameter("kb", part_id = self.IPTG_2_lacI.name, mechanism = mech_b)
            print(self.IPTG_2_lacI.name, kb)

        reactions.append(Reaction([self.IPTG, self.plac_2_lacI], [self.IPTG_2_lacI, self.assembly.dna], k=kb))

        # Unbound transcription plac, kb, ku, ktx
        reactions += mech_tx.update_reactions(dna = self.assembly.dna, component = self, part_id = self.name, \
                                        transcript = self.transcript, protein = self.assembly.protein)

        return reactions

Sort out parameter printer

In [43]:
def param_printer2(x0, parameters):
    inits = pd.DataFrame(list(x0.items()),columns = ['param','value']) 
    

    df = pd.DataFrame(list(parameters.items()),columns = ['param','value']) 
    col_part_id = []
    mech = []
    for val in df['param']:
        if isinstance(val, tuple):
            col_part_id.append(val[0])
            mech.append(val[1])
        else:
            col_part_id.append(' ')
            mech.append(val)

    df['param'] = mech
    df['part'] = col_part_id
    df.reindex(['part','param','value'],axis=1)
    df[''] = ['|']*len(col_part_id)
    
    return pd.concat([df, inits],  axis=1).fillna("")

AND simulation

In [44]:
parameters = {
    ('pttr_phosphate_P_protein_ttrR', 'cooperativity'): 2.0, 
    ('pttr_phosphate_P_protein_ttrR', 'ktx'): 0.17025,  #this is the transcription rate
    ('pttr_phosphate_P_protein_ttrR', 'ku'): 11.01321586, #this is the polymerase unbinding rate

    ('transcription_mm',  'ktx'): .01,    #These are the parameters for transcription leak
    ('transcription_mm','kb'):100, #default binding rate!
    #('ttrS-P',"kdephos"):1000, #auto dephosphorylation
    ('ligand_tt_protein_ttrS-P','kdephos'):1000, #dephosphorylation rate of ttrS
    #('ttrS-P',"kphos"):10, #ttrs phosphorylates itself automatically. This is a proxy for sensing
    ('ligand_tt_protein_ttrS-P','kphos'):10, #equivalent to kphos, from before
    #('ttrS-ttrR','ktransfer'):5, #rate of transfer from ttrS to ttrR
    ('ligand_tt_protein_ttrS-ttrR', 'ktransfer'):5, #this is the new transfer rate
    ('ttrR-P',"kdephos"):1000, #auto dephosphorylation
    "kb":100, "ku":10, "ktx":.05, "ktl":.2, "kdeg":3, 'cooperativity':2, #some default parameters
    
    ('phrpL_hrpR', 'cooperativity'): 1.0, #default is 2, but these proteins bind at cooperativity = 1
    ('phrpL_hrpS', 'cooperativity'): 1.0, 
    ('phrpL_hrpR_hrpS_RNAP', 'ktx'): 0.17025,  #this is the transcription rate
    ('phrpL_hrpR_hrpS_RNAP', 'ku'): 11.01321586, #this is the polymerase unbinding rate

    ('translation_mm', 'B0030', 'ku'): 888, #Unbinding
    ('translation_mm', 'B0030', 'ktl'): .15, #Translation Rate

    ('transcription_mm',  'ktx'): .01,    #These are the parameters for transcription leak
    ('transcription_mm', 'ku'): 100,
    
    ('phrpL_hrpR','ku'):100,
    ('phrpL_hrpS','ku'):100,
    ('phrpL_hrpS_RNAP', 'ku'): 1000.,#we'll say that leak has 100x the unbinding rate
    ('phrpL_hrpR_RNAP', 'ku'): 1000,
    ('phrpL','ku'):100000,
    
    ('binding', 'pLac_lacI','ktx'):0.001,
    ('plac_lacI','kb'):10,
    ('plac_lacI','ku'):0.001,
    
    ('binding', 'plac_hrpS_2x_lacI','kb'):8888,
    ("simple_transcription", 'plac','ktx'):0.1,
    
    #Make a lot of repressor
    ('translation_mm', 'rbsL', 'ku'): 0.0001, #Unbinding
    ('translation_mm', 'rbsL', 'ktl'): 2, #Translation Rate
    ('plac_hrpS','kb'):100,
    ('plac_hrpS','ku'):0.0001
}


#P_rbsS_ttrS
constit_ttrS = DNAassembly(name = "constit_ttrS", promoter = "P", rbs = "rbsS", protein = "ttrS")
#P_rbsR_ttrR
constit_ttrR = DNAassembly(name = "constit_ttrR", promoter = "P", rbs = "rbsR", protein = "ttrR")
#P_arl_lacI
constit_lacI = DNAassembly(name = "constit_lacI", promoter = "P", rbs = "rbsL", protein = "lacI")
#P_arl_lacI
constit_mScarlet = DNAassembly(name = "constit_mScarlet", promoter = "P", rbs = "B0034", protein = "mScarlet")

# plac_hrpS
IPTG = Species("IPTG", material_type='inducer')
#IPTG = Protein("IPTG")

lac_p = LacIPromoter2(name="pLac", IPTG = "IPTG", lacI = "lacI", leak=False)
plac_hrpS = DNAassembly(name = "plac_hrpS", promoter=lac_p, rbs="rbsH",protein="hrpS")

# pttr_hrpR
tt = Species("tt",material_type="ligand")
ttrSbound = ChemicalComplex([Protein("ttrS"),tt])
ttrS3 = PhosphoTransferase(ttrSbound.species,targets=["ttrR"])
ttrR = PhosphoProtein("ttrR")
# Need to change from UTR1
pttr_hrpR = DNAassembly(name = "pttr_hrpR",promoter=RegulatedPromoter("pttr",[ttrR.get_phosphorylated_species()],leak=False),rbs="UTR1",protein="hrpR")

hrpL_gfp = DNAassembly("hrpL_gfp",promoter=CombinatorialPromoter("phrpL",["hrpR","hrpS"], tx_capable_list = [['hrpS', 'hrpR'], ['hrpS'], ['hrpR']], leak=False),rbs="B0030",protein="GFP")

extract_1_TXTL = TxTlExtract(name = "e coli", components = [pttr_hrpR,ttrSbound,ttrS3,ttrR,constit_ttrS,constit_ttrR,
                                                                     constit_lacI, plac_hrpS, hrpL_gfp, constit_mScarlet],parameters=parameters)

CRN_extract_1 = extract_1_TXTL.compile_crn()
r
inducer_IPTG_2x_protein_lacI 100
In [45]:
param_printer2(x0, parameters)
Out[45]:
param value part param value
0 cooperativity 2.000000 pttr_phosphate_P_protein_ttrR | dna_comb_promoter 5
1 ktx 0.170250 pttr_phosphate_P_protein_ttrR | protein_hrpS 1
2 ku 11.013216 pttr_phosphate_P_protein_ttrR | protein_hrpR 1
3 ktx 0.010000 transcription_mm |
4 kb 100.000000 transcription_mm |
5 kdephos 1000.000000 ligand_tt_protein_ttrS-P |
6 kphos 10.000000 ligand_tt_protein_ttrS-P |
7 ktransfer 5.000000 ligand_tt_protein_ttrS-ttrR |
8 kdephos 1000.000000 ttrR-P |
9 kb 100.000000 |
10 ku 10.000000 |
11 ktx 0.050000 |
12 ktl 0.200000 |
13 kdeg 3.000000 |
14 cooperativity 2.000000 |
15 cooperativity 1.000000 phrpL_hrpR |
16 cooperativity 1.000000 phrpL_hrpS |
17 ktx 0.170250 phrpL_hrpR_hrpS_RNAP |
18 ku 11.013216 phrpL_hrpR_hrpS_RNAP |
19 B0030 888.000000 translation_mm |
20 B0030 0.150000 translation_mm |
21 ku 100.000000 transcription_mm |
22 ku 100.000000 phrpL_hrpR |
23 ku 100.000000 phrpL_hrpS |
24 ku 1000.000000 phrpL_hrpS_RNAP |
25 ku 1000.000000 phrpL_hrpR_RNAP |
26 ku 100000.000000 phrpL |
27 pLac_lacI 0.001000 binding |
28 kb 10.000000 plac_lacI |
29 ku 0.001000 plac_lacI |
30 plac_hrpS_2x_lacI 8888.000000 binding |
31 plac 0.100000 simple_transcription |
32 rbsL 0.000100 translation_mm |
33 rbsL 2.000000 translation_mm |
34 kb 100.000000 plac_hrpS |
35 ku 0.000100 plac_hrpS |
In [46]:
x0 = {"dna_pttr_hrpR":1,
      "dna_constit_ttrS":1,
      "dna_constit_ttrR":1,
      "dna_constit_lacI":1,
      "dna_plac_hrpS":1,
      "dna_hrpL_gfp":1,
      "dna_constit_mScarlet":5,
      "phosphate_P":100,
      "protein_RNAP":15,
      "protein_RNAase":30,
      "protein_Ribo":120,
     "protein_lacI": 5}
In [47]:
IPTG_um_list = np.logspace(-3, 2, 6).round(decimals=6)
tt_um_list = np.logspace(-3, 2, 6).round(decimals=6)

fig, heat_map = simulate_heatmap('ligand_tt',tt_um_list,'inducer_IPTG',IPTG_um_list,x0=x0,normalize=True,title='Compiled Full AND gate \n lacI Spike', CRN=CRN_extract_1)

heat_map.set(xlabel='IPTG (uM)', ylabel='Tetrathionate (uM)')
heat_map.invert_yaxis()

plt.yticks(fontsize=18, fontstyle='normal', rotation = 0)
fig.savefig("./figures/fig6/fig6d.png")
plt.show()
/Users/lianamerk/anaconda3/envs/python3.7-bioscrape/lib/python3.7/site-packages/biocrnpyler-0.1-py3.7.egg/biocrnpyler/chemical_reaction_network.py:870: UserWarning: The following species are uninitialized and their value has been defaulted to 0: inducer_IPTG, 
/Users/lianamerk/anaconda3/envs/python3.7-bioscrape/lib/python3.7/site-packages/scipy/integrate/odepack.py:248: ODEintWarning: Excess work done on this call (perhaps wrong Dfun type). Run with full_output = 1 to get quantitative information.
  warnings.warn(warning_msg, ODEintWarning)
odeint failed with mxstep=500...odeint failed with mxstep=500...odeint failed with mxstep=500...odeint failed with mxstep=500...odeint failed with mxstep=500...odeint failed with mxstep=500...
In [48]:
x0 = {"dna_pttr_hrpR":1,
      "dna_constit_ttrS":1,
      "dna_constit_ttrR":1,
      "dna_constit_lacI":1,
      "dna_plac_hrpS":1,
      "dna_hrpL_gfp":1,
      "dna_constit_mScarlet":5,
      "phosphate_P":100,
      "protein_RNAP":15,
      "protein_RNAase":30,
      "protein_Ribo":120,
     "protein_lacI": 0}
In [49]:
IPTG_um_list = np.logspace(-3, 2, 6).round(decimals=6)
tt_um_list = np.logspace(-3, 2, 6).round(decimals=6)

fig, heat_map = simulate_heatmap('ligand_tt',tt_um_list,'inducer_IPTG',IPTG_um_list,x0=x0,normalize=True,title='Compiled Full AND gate \n No lacI Spike', CRN=CRN_extract_1)

heat_map.set(xlabel='IPTG (uM)', ylabel='Tetrathionate (uM)')
heat_map.invert_yaxis()
plt.yticks(fontsize=18, fontstyle='normal', rotation = 0)
fig.savefig("./figures/fig6/fig6b.png")
plt.show()

Early Circuit Dynamics in no lacI Spike Model

In [50]:
timepoints = np.linspace(0, 200, 1000)

x0 = {"dna_pttr_hrpR":5,
      "dna_constit_ttrS":1,
      "dna_constit_ttrR":1,
      # Cole origin
      "dna_constit_lacI":40,
      "dna_plac_hrpS":1,
      "dna_hrpL_gfp":1,
      "dna_constit_mScarlet":5,
      "phosphate_P":100,
      "protein_RNAP":10,
      "protein_Ribo":90.,
      "ligand_tt": 1,
      "inducer_IPTG": 0.,
      "protein_lacI": 0,}
Re1 = CRN_extract_1.simulate_with_bioscrape(timepoints, initial_condition_dict = x0)

results_melted = pd.melt(Re1, var_name='species', id_vars=['time'])
results_melted.head()
Out[50]:
time species value
0 0.000000 complex_2x_phosphate_P_protein_ttrR_dna_pttr_h... 0.000000e+00
1 0.200200 complex_2x_phosphate_P_protein_ttrR_dna_pttr_h... 5.330060e-13
2 0.400400 complex_2x_phosphate_P_protein_ttrR_dna_pttr_h... 2.037818e-11
3 0.600601 complex_2x_phosphate_P_protein_ttrR_dna_pttr_h... 1.643681e-10
4 0.800801 complex_2x_phosphate_P_protein_ttrR_dna_pttr_h... 7.002289e-10
In [51]:
Re1 = CRN_extract_1.simulate_with_bioscrape(timepoints, initial_condition_dict = x0)
/Users/lianamerk/anaconda3/envs/python3.7-bioscrape/lib/python3.7/site-packages/biocrnpyler-0.1-py3.7.egg/biocrnpyler/chemical_reaction_network.py:870: UserWarning: The following species are uninitialized and their value has been defaulted to 0: inducer_IPTG, 
In [52]:
p = bokeh.plotting.figure(title = 'Early Circuit Dynamics, pLac', x_axis_label = 'Time (min)', y_axis_label = 'Simulated Concentration (uM)', height = 2000, width = 4300)

p.line(timepoints, Re1['complex_dna_plac_hrpS_protein_RNAP'], legend_label = 'Active pLac', color = 'Green', line_width = 14)
p.line(timepoints, Re1['rna_constit_lacI'], legend_label = 'LacI RNA', color = 'orange', line_width = 14)

p.line(timepoints, Re1['protein_hrpS'], legend_label = 'hrpS protein', color='lightsteelblue', line_width = 14)

p.line(timepoints, Re1['complex_dna_hrpL_gfp_protein_hrpR_protein_hrpS_protein_RNAP'], legend_label = 'hrpS, hrpR both bound to pHRPL', color='darkorchid', line_width = 14)
p.line(timepoints, Re1['rna_hrpL_gfp'], legend_label = 'GFP RNA', color='navy', line_width = 14)

p.add_layout(p.legend[0], 'right')
p.legend.title = 'Species'


p.legend.glyph_width = 100
p.legend.padding = 100
p.legend.title_text_font_size = '80pt'
p.legend.label_standoff = 50
p.legend.glyph_height = 100

p.xaxis.major_tick_line_width = 10
p.yaxis.major_tick_line_width = 10

p.xaxis.minor_tick_line_width = 10
p.yaxis.minor_tick_line_width = 10
p.axis.minor_tick_out = 20
p.axis.minor_tick_out = 20

p.axis.major_tick_out = 40
p.axis.major_tick_out = 40
p.yaxis.axis_label_standoff = 50
p.xaxis.axis_label_standoff = 50

p.xaxis.axis_line_width = 8
p.yaxis.axis_line_width = 8

p.yaxis.major_label_standoff = 30
p.xaxis.major_label_standoff = 30

p.ygrid.grid_line_width = 10
p.xgrid.grid_line_width = 10

p.axis.axis_label_text_font_size = "80pt"
p.legend.label_text_font_size = "80pt"
p.title.text_font_size = '100pt'
p.yaxis.major_label_text_font_size = "80pt"
p.xaxis.major_label_text_font_size = "80pt"

#p.y_range = Range1d(0, 0.9)
bokeh.io.show(p)

Computing Environment

In [53]:
%load_ext watermark
%watermark -v -p numpy,pandas,scipy,bioscrape,biocrnpyler,bokeh,jupyterlab,holoviews
CPython 3.7.6
IPython 7.13.0

numpy 1.18.1
pandas 1.0.3
scipy 1.4.1
bioscrape 1.0.0
biocrnpyler unknown
bokeh 2.0.1
jupyterlab 2.1.0
holoviews 1.13.2
/BrvoWdNsxL8+bNWbBgARcvXqR27doAugUdSvcXNQdZnSiEECXUQ62mwI+5aDQabt68yYMHD4CcLdyaNm3KuHHj+OWXX4iLi2P69On06NHD4DSGZ0mSmBBClFDZaAv8mMv169dp3bo1u3fvBkClUhEaGkqNGjUYMGAA/v7+tGrVipkzZ5rtnkrIdKIQQpRQGjMmqSc9eQhs7ruMj6tUqZLeeXhFQZKYEEKUUA+12UUdQpErVtOJ06ZNw8/Pz6A8KSmJ0aNH4+npiaenJxMmTOD27dtPXc+YwrQVQoiikK3g809XbEZikZGRREZG4unpqVeemprKgAEDUKvVDB06FI1Gw+rVq0lISCAyMlL3sp3SesYUpq0QQhQVtRxCUvRJTKPRsGLFCqMv4wGsXbuW5ORkdu7cSZ06dYCcVTGDBg0iKiqKd99916R6hbmHEEIUJy/CSKsgRTqdmJWVhY+PD0uXLs1zWWZ0dDSenp665ALQqlUrnJ2diY6ONrmeMYVpK4QQReWhVlXg559O0UhMo9Hw448/EhcXx+nTp7lx4wapqamULl2acuXK4ezsjJubG61atTLpJbesrCzS09MJDg6mW7duBhtP3r17l6SkJLp06WLQ1s3NjUOHDplUz5jCtBVCiKKk4Z+fpAqSbxK7ffs2a9euZfv27fz9998APHkQ9N27d7ly5QpHjhxhxYoVVK1alR49ejBw4EDKly+f781tbGyIiYmhVCnjYdy4cQPA6AjN3t6e9PR07t27p7iera3tU9/DWFshhChKD7XFam1ekcgzia1atYovv/yS+/fvo9VqKVeuHPXr16du3brY29vrjp9OT0/n1q1bnD9/noSEBK5fv87KlSsJDw9n8ODBjBo1Ks+bW1hYYGGR97+E+/fvA/Dyyy8bXHvppZcAyMjIUFzPWCIqTFshhChKMhLLJ4ktXLgQBwcH+vbtS4cOHXB3d0elKvgXlpiYyL59+4iKimLp0qX5JrGCZGcX/NjSwsJCcb3C3EMIIYqbh1rLog6hyOWZxL788kveeOMNk/+A161bl7p16zJy5EhiY2MLFVzuaM/Y2TS5ZWXLllVcrzD3EEKI4kZGYvkksbZt2xa689dff71Q7atVqwbAzZs3Da6lpKRgZ2dHmTJlFNcrzD2EEKK40cgzsbyTWGBgoMmdqVQqZs2aVaiAHmdnZ4eTkxNnzpwxuHb27FkaNGhgUr3C3EMIIYqbh8h0Yp5JLDIyUu8ZmFarzfeZWO51cyYxAC8vL9avX8+FCxd073EdPXqUS5cuMWTIEJPrFeYeQghRnMhITOF7Yo6Ojs/1pM7HDRs2jO+++46BAwcyePBgsrKyCAsLw83NjR49ephcLykpiVOnTtG0aVNq1KhhUlshhChOZCRWQBLLHV1dv34dAA8PD5o2bYqHhwf16tV7LgFWrFiR8PBwgoKCWLJkCdbW1nTq1IkJEybo7WmotF58fDyTJ08mKChIl8SUthVCiOJERmL5JLGgoCDi4+M5ceIEV65c4a+//uL69evs2rULAFtbW5o0aUKTJk3w8PCgUaNGuveqntbBgweNlteuXZtVq1YV2F5JPV9fX3x9fZ/6HkIIUVzIEvt8kpiPjw8+Pj5Azsq9+Ph43efChQukpaVx+PBhjhw5ktNRqVK4ubmxefPm5xO5EEK84DTF6zStIqHomZi9vT3dunWjW7duQM7RJREREaxbt467d+8C8PDhQ37++ednF6kQQgg9D7VFfhBJkVP8Gzh37hxxcXG6Kcbc5PU4eSlYCCGeH80LsEt9QfJMYr/++ivx8fHExcVx6tQp7t27B+hvAOzo6EjTpk11HxcXl2cfsRBCCEBGYpBPEuvVq5fuvTCtVoulpSX16tXTS1qOjo7PLVAhhBD6smXbKWXTidWqVaNx48aULVsWtVrNsWPHOHbsmEG9Z/GysxBCCONK2hL7v//+m+3bt/PXX39Rs2ZNunfvTsWKFQvVp0r75AFh/8/V1RWVSlXgTh25cuudO3euUAEJIYRQJvT3DgXWGe1q/NWl5+38+fP069dPbz1FuXLlWLZsmUmHKT8pz5FYkyZNFCUvoc/r9c+KOgSTxMQGlqiYY2ID6fyv2UUdhkn2/XcabbrPL+owFDuyYzwAf14tOY8LajldJzv5+WzAYC4WVf8odB8laYn94sWLKVu2LMuWLaNhw4ZcvHiRqVOnMmvWLHbs2PHU/eaZxCIiIp66UyGEEM+euV92zs7OJjQ0lMjISNLS0vDw8GDGjBnUrFnTaP2UlBSCgoI4evQoAC1btmTy5MlGtyk8ffo0EyZM0I26XnvtNaZOnUq/fv24ffv2U08rlpw0LoQQQk+21qLAjymWLVtGREQEs2fP5ptvvsHS0pIhQ4YYPW8RYMyYMVy/fp01a9bw9ddfk5yczMiRI43WvXPnDtWrV9crc3FxQavVcuvWLZPifJzi9Znvv/9+gXVUKhXh4eFPHYwQQgjlzDkSU6vVrFmzhoCAAN15ksHBwbRu3Zo9e/bg7e2tV//27ducPn2aFStW4ObmBsDw4cMZNWoUf//9N5UqVdKrr9FoDA5Zfvnll3O+x8OHTx234iR28uRJ3UKPxz2+DF+eoQkhxPNjztWJ586dIyMjg5YtW+rKbGxsqF+/PidOnDBIYmXKlKFMmTJERUXh6emJSqVi165d1KpVi/Lly5stroIoTmJPLvTQarWo1Wpu3LjBrVu3eO2112jatOkzCVIIIYQhJSOxtLQ00tLSDMrt7Oyws7PT/Xzjxg0AqlSpolfPwcFBd5LJ46ytrQkKCmLmzJk0a9YMlUpF5cqVCQ8Px9LSeFzZ2dlkZ2frftZoNEbLAYNRW14UJ7H8FnqEhoaycuVKpkyZorQ7IYQQhZStYNupdevWERoaalA+evRo/P39dT9nZmYCGBw/ZWVlhVqtNmiv1Wo5e/Ys7u7uDB8+HI1GQ0hICKNGjWLz5s3Y2toatPHz8zMa47vvvqv3s0ql4uzZswV+NzAhieXngw8+YNmyZYSEhLBx40ZzdCmEEKIASkZiAwYM0J1I8rjHR2GQM7KCnGdjjycytVpNmTJlDNrv3r2bjRs3cujQIV3CWrFiBe3bt2fLli0MGTJEr/7o0aML/kJPwSxJ7KeffkKr1XLmzBlzdFegadOmcfnyZTZs2KBXnpSUxLx584iLiwOgXbt2TJo0SdHSzcK0FUKIopCtYIH5k9OGecndRjAlJQUbGxtdeUpKCnXr1jWof/LkSWrWrKk34ipXrhzOzs5cvnzZoH6RJzFjqxOzs7O5f/8+Fy9eRKVSGaxGeRYiIyOJjIzE09NTrzw1NZUBAwagVqsZOnQoGo2G1atXk5CQQGRkZL4nNBemrRBCFBVz7mLv6uqKjY0NcXFx1K5dG4D09HTOnj1L3759DepXrVqVK1eukJmZqVtlmJGRwdWrV3nrrbcM6i9atIg+ffqYfc/dQq9OfNzAgQPNEZNRGo2GFStWGJ3bBVi7di3Jycns3LmTOnXqAODu7s6gQYOIiooymHM1V1shhCgqj7LNt8TeysqKfv36ERwcTOXKlXFycmLhwoVUqVIFLy8vNBoNt2/fxtbWFmtra7y9vVm9ejVjx47l448/BiAkJITSpUvTs2dPg/5XrVpFu3btdElMq9USEBDA2LFjcXJyeuq4n3p1Yq7SpUtjb2+Pl5cXXl5eTx1IfrKysujVqxcJCQl4e3sTGxtrUCc6OhpPT09dEgJo1aoVzs7OREdH55uICtNWCCGKisbMu9iPGTMGjUbD9OnTyczMxMPDg7CwMKysrLh69SodO3YkKCgIX19fHBwc2LRpE/Pnz9cNYDw8PIiIiKBcuVzsm5YAACAASURBVHIGfT85AMrOziY6OprBgwc/nyRWlNtQZWVlkZ6eTnBwMN26daNDB/1NL+/evUtSUhJdunQxaOvm5sahQ4fy7LswbYUQoiiZcyQGYGlpSUBAAAEBAQbXnJycSEhI0CurU6cOK1euNGsMpsrzqaCxB3OmSk5OLnQfkPPCXUxMDN26dTN6Pa/3GwDs7e1JT0/XHeppzrZCCFGUslEV+PmnyzOJdenSBT8/P7799ltSU1MVd5ienk5MTAyjRo2iU6dO5gnSwoJSpfIeNN6/fx/43xYmj3vppZeAnAeO5m4rhBBF6WG2ZYGff7o8M0OPHj3YsWMHJ06cwMLCgoYNG9KgQQPq1q2Lvb09ZcuWxdLSkvv375OcnMyFCxc4efIkCQkJujev8xo5mduTb3obk9fb34VpK4QQRUnJy87FyeM7c+S3WweYYceOefPm4efnx8qVK/nPf/7DTz/9xM8//5xvZ1qtllKlStGpUydGjhxJ/fr1FQVRWGXLlgUwutNyblluHXO2FUKIovSohJ3sbGzHDmML58y2Y0eDBg0IDQ3l+vXr7N+/n6NHj3L69Gnu3LmjV8/BwYHGjRvTokULunbt+txfEK5WrRoAN2/eNLiWkpKCnZ2d0TfOC9tWCCGKkqlHrRSlIn3Z2dHRET8/P10WffToEampqWi1WipUqEDp0qWfSXBK2dnZ4eTkZHTHkLNnz9KgQYNn0lYIIYpSSRqJPask9lS/gVKlSmFvb4+Dg0ORJ7BcXl5exMbGcuHCBV3Z0aNHuXTpUoHP5grTVgghikq2VlXgp7j766+/9D4PHjwwqb1Z9k4sDoYNG8Z3333HwIEDGTx4MFlZWYSFheHm5kaPHj109ZKSkjh16hRNmzalRo0aJrUVQoji5FF2yRmJQc7pzosXL6ZChQq6F6s7dOigt5FGhw4dWLZsmeI+S9ZvIB8VK1YkPDwcV1dXlixZwrp16+jUqZPubfNc8fHxTJgwgfj4eJPbCiFEcVKSRmL37t3jvffeIyoqSvf6Uq4RI0Ywa9Ys3n//fQ4ePMjp06cV91siR2IHDx40Wl67dm1WrVqVb1tfX198fX2fqq0QQhQnJell5q+//po7d+7w7bffUqtWLb1rnTp1ws3NDYDjx48TGRlJkyZNFPX7jxmJCSHEi+ZRtkWBn+Ji3759+Pn5GSSwJ/Xo0YMTJ04o7vepv6FGozFYai+EEOL5KUnTiVevXjUYXalUKl5++WUsLf+3s8hrr72m2w5QCZOSWFZWFl9++SXdu3fH3d2dVq1aATnzmb/99pspXQkhhCgkTbZFgZ/i4vFElcvCwoLTp0/j6uqqK3v48KHulGklFD8Tu3v3Lv369SMxMVG3pb5KpeLGjRscOnSI+Ph4wsPDee211xTfXAghxNMrSc/EqlWrxrlz5/jXv/6Vb71ff/1Vt3JcCcVpevHixZw/f56KFSvqHX6ZmZmJg4MDGRkZLF26VPGNhRBCFE5JGom1b9+eiIgIMjMz86xz7949IiMjTTqbUvE33LdvHyqVirlz5zJo0CBdea1atQgKCkKr1Zq0LFIIIUThlKRnYv369SMjI4P+/fsbnEsGOe/wDh8+HJVKRe/evRX3q3g6MXcRR/Xq1Q2uOTg4AJj8pvU/UUxsYFGHYLKSFvO+/04r6hBMdmTH+KIOwWS1nK4XdQgmsaj6R1GH8NwVp5FWQezt7Vm6dCkff/wx3t7evPrqqzg7O6NSqbh69Spnz56lYsWKhIaGYmdnp7hfxUmsVq1aJCYmsnbtWj744ANd+aNHj/jyyy91dV50Xq9/VtQhmCQmNrBExRwTG0jnVrOLOgyT7Ds6jbZvzy/qMBQ7vDMn4Z668koRR6Jc01eu8PdfT3/EfVGoVO1qofv4/+UJJUazZs3YtWsXERER7Nu3j6NHj6LRaKhevTrDhw/Hz8+PSpUqmdSn4iTWr18/ZsyYQWRkJNu3b9dtE9KsWTOysrJQqVRGt9QXQgjxbGhK0AbAucqXL8/IkSMZOXKkWfpT/Bt47733GDZsGJAz+tJqtWi1Wh48eICFhQX9+/enT58+ZglKCCFEwUrSM7FnxaRtp8aNG0fPnj3Zv38/V65cwdLSEicnJ7y8vExaEimEEKLwStp04rOgOInt3LkTyDm2ZOjQoXrXkpOTWbt2LeXKlcPHx8e8EQohhDAquwQt7HhWFCex8ePHY2FhQYsWLXSrEXPdv3+fuXPnYm9vL0lMCCGekxdhurAgeSaxxMREPv/8c70yrVbLJ598onc8iVar5a+//gIgPT39GYUphBDiSdnZ5k1i2dnZhIaGEhkZSVpaGh4eHsyYMYOaNWsarf/w4UOWLFlCVFQU9+7do0GDBkydOvW57tyUZxKrW7culpaW/PDDD6hUKt1qxJMnTxrUzd2GqmHDhs8ozKczbdo0Ll++zIYNGwqsm5SUxLx584iLiwOgXbt2TJo0iYoVKz7rMIUQ4qlozTwSW7ZsGREREcydO5cqVaqwcOFChgwZQnR0tMEZYAAzZ87kwIEDzJ07lxo1arB48WKGDh3Knj17DN71io2NNSmW119/XVG9fKcTAwMD+eSTT9BqtZw5cwaVSkW9evUoVep/zVQqFVZWVtStW1e3erE4iIyMJDIyEk9PzwLrpqamMmDAANRqNUOHDkWj0bB69WoSEhKIjIyUgzGFEMWSOacT1Wo1a9asISAggLZt2wIQHBxM69at2bNnD97e3nr1k5KS2Lp1K8uWLaNdu3YAzJkzhx49evDLL7/QunVrvfqDBg1CpVLp7b2by1jZuXPnFMWdbxJ75ZVX2Lp1KwBt2rRBpVKxevVqKleurKjzoqDRaFixYgWhoaGK26xdu5bk5GR27txJnTp1AHB3d2fQoEFERUXJ+29CiOLJjKsTz507R0ZGBi1bttSV2djYUL9+fU6cOGGQxH788UfKli1L+/btdWW2trZ5Hlq8fv163T/fuHGDqVOn0r17d/7973/j4OBAamoq+/fvZ8uWLQaPsvKjeGHHkSNHCqyTnJxM1apVFd/c3LKysujVqxcJCQl4e3srHr5GR0fj6empS2AArVq1wtnZmejoaEliQohiSckzsbS0NNLS0gzK7ezs9Kb8cs/wqlKlil49BwcHrl833ILszz//xMnJiUOHDrFixQquX79O/fr1mTRpkt7f0lyPz4oNGTKE3r17M2XKFL06Hh4elC5dmq+//pouXboU+N3AxPfEoqKi+OGHH7h79y7Z2dm68kePHpGamsrFixc5c+aMKV2aVVZWFunp6QQHB9OtWzc6dOhQYJu7d++SlJRk9Bfm5ubGoUOHnkGkQghReEqeia1bt87ozNTo0aPx9/fX/Zy7u/yTj0+srKxQq9UG7dPT07l27RohISGMHz+e8uXLs3LlSvr27Ut0dHS+M3YnTpxgwIABRq81a9aMdevWFfi9cilOYps3b2bmzJm6Oc0n5zZzy4qSjY0NMTExes/sCpLXf31AzoaV6enp3Lt3D1tbW7PFKYQQ5qBVMBIbMGCA0Vefnlx4kXsQpVqt1ktkarWaMmXKGLQvXbo06enpLFiwABcXFwAWLVpE27Zt2bZtm94eu0+yt7fn5MmTtGnTxuDakSNHcHR0LPB75VL8137Lli0AODo64uDgwM8//0y7du149OgRsbGxZGdnM3fuXMU3fhYsLCywsDDt5b/79+8D8PLLLxtcy12Nk5GRIUlMCFH8KHgm9uS0YV5yE0dKSgo2Nja68pSUFOrWrWtQv2rVqqhUKl599VVdmbW1NTVq1ODq1fw3N37vvfcIDg7mwYMHdOrUiYoVK3Lz5k12795NZGQkn376acFf7P8pTmKXL19GpVIRHBxMuXLlePPNN+nZsyedO3dm9erVzJ8/n9jYWHr06KH45sXB49OieTE1MQohxPNgziX2rq6u2NjYEBcXR+3atYGcKcOzZ8/St29fg/rNmjVDq9Xy22+/0ahRIyDnOK68Hs88btiwYaSlpbF27Vrdgg+tVsvLL7/M+PHjTVqHoDiJaTQaACpVqoSTkxMVKlTg9OnTdO7cmY4dOzJ//nyOHz+u+MbFRdmyZYGc52lPyi3LrSOEEMWJkulEpaysrOjXrx/BwcFUrlwZJycnFi5cSJUqVfDy8kKj0XD79m1sbW2xtramWbNmtGrViokTJzJr1iwqVKjAkiVLUKlU+Pr6Fni/cePG8cEHH/DTTz9x584dKlasSJMmTYzOiuVH8RAj95lRREQEGo2GRo0asXfvXm7cuEFMTAzwv4MzS5Jq1aoBcPPmTYNrKSkp2NnZGZ0PFkKIIqdV8DHBmDFj6NWrF9OnT6dPnz5otVrCwsKwsrLi+vXrtG7dmt27d+vqh4aG0rJlS/z9/enZsydpaWmsX79e8ZlgNjY21KlTBycnJxo3bvxU6yoUj8S6du3KypUrWbNmDe+88w5vvPEGhw8f1r3kplKpcHZ2NjmAomZnZ4eTk5PRVZVnz56lQYMGRRCVEEIoYOYdOywtLQkICCAgIMDgmpOTEwkJCXplZcuWZcaMGcyYMcPkex06dIh58+bx559/olKpiIyMZOnSpTg4ODBz5kzFj3EUj8RGjRpF9+7dKVOmDM7Ozvj6+lKzZk3duWKWlpZ89NFHJn+R4sDLy4vY2FguXLigKzt69CiXLl2iW7duRRiZEELkw8wjsefl8OHDjBo1imrVqjF9+nTd2oQWLVqwbds2wsLCFPeleCRmZWXFF198we3btwEoU6YM27ZtY9euXWRkZNCmTRujK1iKm6SkJE6dOkXTpk11Z6ANGzaM7777joEDBzJ48GCysrIICwvDzc2txC1UEUK8OMz5TOx5WrJkCV26dCE4OBiNRqNbjTho0CDu3LnD9u3bGT58uKK+TF529/iGuDY2NvTu3ZvBgweXiAQGEB8fz4QJE4iPj9eVVaxYkfDwcFxdXVmyZAnr1q2jU6dOurlgIYQolkroSOz8+fN5DhBatmxpdIeQvJi0Y8fZs2fZvXs3t2/f5tGjRwbXVSoV8+bNM6XLZ8rYHl6+vr5GV87Url2bVatWPY+whBDCLFQldCRmZ2enO8LrSVevXjXpvVzFSWzHjh1MnDgxz+u5O3YUpyQmhBD/aMV0pFWQjh07snTpUurVq0eTJk2AnEHQtWvX+OqrrxRtGZhLcRJbuXKlbpupl156CVtbWywtLU0MXQghhNmU0JHYuHHj+OWXX/Dz86NChQoAfPzxxyQnJ1OjRg0++eQTxX0pTmJXr15FpVIxdOhQPvnkkyLfJ1EIIV54JXQkZmdnxzfffMN3333HsWPHSE1NxdbWlgEDBuDr62vSC8+Kk1itWrU4f/48Pj4+ksCEEKI4KKEjMchZ8d6rVy969eqlV56enk58fDzNmzdX1I/i1YljxowBco5jEUIIUfRU2oI/xdFrr71GYGCg0QWCCQkJ9O/fX3Ffikdily5dolGjRqxatYqDBw9Su3Zt3db9uWRhhxBCPEfFNEkVRKvVsnXrVi5dusSSJUv0Xt0ylUqbu1qjAK6urvmeG5Zbfu7cuacORgghhHK1lywssM7FMeOeQySmcXV1Zfr06bqdOZYvX46rqysAJ0+epF+/fopzieKRmIODgzwLU6Cry6SiDsEkexLmlqiYS1q8kBNzl+bKz0cqanvjc/bBq714URFHotzFjz7BeWnBf9CLk0v+ZkguZt478XlydXVly5YtjBw5kj59+jB37ly6dOlicp5RnMSOHDlicpBCCCGeoYKPQyzWKleuTHh4OBMmTODjjz/mww8/1G0qr5RJO3bkyszM5PLlywA4OzvrTkAWQgjx/BTXhRumeOmll1i8eDHBwcGEhoby448/mtTepCR27949Pv/8c3bt2qVbVVK6dGm8vb2ZMGGC3pHWQgghnrESPhJ73NixY3F2diYwMNCkdoqTWHp6On379iUxMZHH14Ko1WoiIyP5+eef2bRpk5yCLIQQz0lJHYkdOHAAe3t7g3Jvb29q1qzJf//7X8V9KU5iK1as4Pz581haWuLt7U2LFi3Izs4mPj6eqKgo/vjjD7788kuTtgsRQghRCCXoZefs7GzdQZeOjo66sie5u7vj7u6uuF/FSWzv3r2oVCrGjh3L0KFDdeXe3t44OzuzYMECdu/eLUlMCCGek5I0EnNzc2PTpk00adKE+vXr57sKUaVScfbsWUX9Kk5iN27cAKB9+/YG19q3b8+CBQtISUlR2t0zk5SUxLx584iLiwOgXbt2TJo0yaSX6aZNm8bly5fZsGHDswpTCCEKTVWCnol9+OGHuhHYhx9+aLZXthQnsUqVKnHjxg1++eUX6tSpo3ft119/BXKWSxal1NRUBgwYgFqtZujQoWg0GlavXk1CQgKRkZGKDriMjIwkMjIST0/P5xCxEEIUQgkaiY0ePVr3z/7+/mbrV3ESa9OmDVu2bGH27NncuXNH90c+Li6OZcuWoVKpeOONN8wW2NNYu3YtycnJ7Ny5U5do3d3dGTRoEFFRUbz77rt5ttVoNKxYsYLQ0NDnFa4QQhROCUpiSUlJJtWvUaOGonqKk9ioUaPYu3cvd+/e5YsvvtC7ptVqsbOzY8SIESYFaW7R0dF4enrqjRRbtWqFs7Mz0dHReSaxrKwsevXqRUJCAt7e3sTGxj6vkIUQ4qmVpGdinTt3NmkK0ezbTlWtWpWNGzcyZcoUfvnlF71r9erVY968ebr5zqJw9+5dkpKS6NKli8E1Nzc3Dh06lGfbrKws0tPTCQ4Oplu3biadKiqEEEWmBCWxoKCgZ9KvSS87161bly1btnD+/HndcvtatWrh4uLyTIIzRe7CkypVqhhcs7e3Jz09nXv37mFra2tw3cbGhpiYGEqVeqoNTIQQokiYe2FHdnY2oaGhREZGkpaWhoeHBzNmzKBmzZoFtt25cycBAQHExMQYre/j42PeYP/fU/3VfvXVV3n11VfNHUuh3L9/H8DoiaC522JlZGQYTWIWFha69xeEEKLEMPNIbNmyZURERDB37lyqVKnCwoULGTJkCNHR0fluL3jt2jU+/dS0Ta4zMzM5d+4carVat4GGVqslIyODkydPMnHiREX95JnEJkyYgEqlYtq0adja2jJhwoQCOyvK88SMvTT3JElUQoh/EnOOxNRqNWvWrCEgIIC2bdsCEBwcTOvWrdmzZw/e3t5G22VnZzN+/Hjc3Nw4duyYonvFxsby8ccfk5aWZvS6jY1N4ZPYjh07UKlUfPLJJ9ja2up+zkvueWJFlcRyt7vKysoyuJZbJltiCSH+Scy5sOPcuXNkZGTQsmVLXZmNjQ3169fnxIkTeSaxlStX8vDhQ0aPHq04iS1evJjy5cvz2WefsWvXLlQqFb6+vhw+fJjNmzfz1VdfKY47zySWe35Y7nOi4n6eWLVq1QC4efOmwbWUlBTs7OwoU6bM8w5LCCGeHTOOxPJaV+Dg4MD169eNtvnll19Ys2YNW7du1bVXIiEhgU8//RQvLy8yMjKIiIigbdu2tG3blgcPHrB8+XJWrVqlqK88k9iT54cV9/PE7OzscHJy4syZMwbXzp49S4MGDYogKiGEeHaUjMTS0tKMTtvZ2dlhZ2en+zkzMxPAYFMIKysr1Gq1QfuMjAwCAgIICAigVq1aJiWx7OxsHBwcgJzjvP744w/dtS5dujB16lTFff2jHhJ5eXkRGxvLhQsXdGVHjx7l0qVLdOvWrQgjE0KIZyC74M+6devo2LGjwWfdunV6XVlbWwMYJCy1Wm10Fmv27NnUqlWL3r17mxz2K6+8QkJCAgC1atUiMzOTS5cuAfDo0SPdQj0lTFqdmJiYSHR0NCkpKTx69EjvSBYo2oUdAMOGDeO7775j4MCBDB48mKysLMLCwnBzc6NHjx5Azlvjp06domnTporfCBdCiOJIyUhswIABRpe3Pz4Kg//tLJ+SkqJ3NmRKSgp169Y1aL9t2zasrKxo0qQJkLPrEUCPHj3o3r07s2bNyjOmt99+m+DgYLKzsxk0aBCNGzdm1qxZ9O3bl5UrVxq9X14UJ7GYmBjGjh2b5yrAol7YAVCxYkXCw8MJCgpiyZIlWFtb06lTJyZMmKAbIsfHxzN58mSCgoIkiQkhSjYFSezJacO8uLq6YmNjQ1xcHLVr1wZyzpE8e/Ysffv2NagfExOj9/PPP//M+PHjWbFiBfXq1cv3XkOHDuXOnTv89ttvAAQGBjJkyBD8/f2xtbVl+fLlBX+x/6c4iS1dulSXacuUKUPZsmWL5ZL12rVr5/tA0NfXF19f33z7OHjwoLnDEkIIszPnEnsrKyv69etHcHAwlStXxsnJiYULF1KlShW8vLzQaDTcvn0bW1tbrK2tDV5oTk5OBnIW2VWqVCnfe1lYWOi9tuXm5sb+/fu5ePEitWvX1hsJFkRxErt8+TIqlYoPP/xQbzdiIYQQRcTMLzuPGTMGjUbD9OnTyczMxMPDg7CwMKysrLh69SodO3YkKCiowIGAKXJn98qUKaNbgPf4AZoFUZzEcleQvPXWW08RphBCCHMz9wbAlpaWuhWHT3JyctItxjCmRYsW+V5/3F9//cWMGTM4deoUGRkZBtefyaGYn3zyCaNGjWLHjh189NFHSpsJIYR4RkrSLvaPCwwM5NSpU/j4+FChQoVCvYOc77ZTT3J0dGTlypXExMRQp04d3ZLMXEW9sEMIIV4oJehk58edPn2aKVOm0KtXr0L3VeC2U0/SarVcvHiRixcvGpRLEhNCiOenpI7EypcvT8WKFc3SV4HbTgkhhCiezH0Uy/PSv39/wsLCaN68uaLl//lRaZ98Y1kIIUSJ0GRUcIF1Ti8f+xwiMU1GRgbvvvsu165dw9nZ2eAILZVKRXh4uKK+FC/sWLlyJSqVCj8/P4MtSK5cucKGDRuwt7dn+PDhSrv8R+pi/X5Rh2CSvQ82lqiY9z7YSJeX/Yo6DJPszdxAF5sBRR2GYnvTc7Yj6vpqwccvFRd7zn9B13rKju4oLvb8UfhHLyV1JBYYGEhiYiK1atUq9OkiipNYSEgIKpUKHx8fgyT24MEDNmzYQIUKFV74JCaEEM9NCZ1HO3jwIGPHjuWDDz4odF95JrHExEQ+/PBDg/K+fftiaWmp+1mr1ZKamgrkbNwohBDi+SipCzusra1xc3MzS195JrG6devSqFEjdu7cCaBb5HH16tU8O/Pw8DBLUEIIIQqmyi6ZWczHx4eNGzfSsmVL3ZmVTyvf1hMnTkStVqPVaomJiUGlUtG2bVteeuklXR2VSoWVlRWvvvoq7777bqGCEUIIYYKSmcOwtrYmPj6e9u3b07BhQ8qWLau3Gt6U17XyTWKVK1dm8eLFAPTp0weVSkVQUJDZ1vcLIYR4eiV1YUdUVJRuaf3vv/9ucN2U17sUj+MiIiIKrHPy5EmZUhRCiOekpD4Ti4iIoEqVKmbpS3ESy87OZvny5fz444/cvXtX71yxR48ecefOHTIzMxVv2miqpKQk5s2bR1xcHADt2rVj0qRJBY4Kn6bdtGnTuHz5Mhs2bDDfFxBCCDMrqSOxnj17EhAQgLe3d6H7UpzEwsLCCA0NRaVSGT3ROXfbqWchNTWVAQMGoFarGTp0KBqNhtWrV5OQkEBkZKTuwEtztIuMjCQyMhJPT89n8l2EEMJsSuhI7NGjR89+26kn5e6l+Nprr1GtWjUOHDhAr169UKlU7NixgwcPHrBy5UqzBPWktWvXkpyczM6dO6lTpw4A7u7uDBo0iKioqDwXlJjSTqPRsGLFCkJDQ5/JdxBCCHMrqasTBw8ezBdffIGFhQWvvvoq9vb2BnXMfp7YtWvXAPj8888pV64c+/fvp0WLFrz11ls0bNiQadOmsXXrVtq2bau0S8Wio6Px9PTUJSKAVq1a4ezsTHR0dJ5JTGm7rKwsevXqRUJCAt7e3sTGxpr9OwghhLmV1GdiW7du5a+//mLYsGFGrz+T88RypwrLlCmDo6MjDg4OnD59mrfeeovmzZsD8NNPPyntTrG7d++SlJREly5dDK65ublx6NChQrfLysoiPT2d4OBgunXrRocOHcwVvhBCPDMqTVFH8HS6d+9utr4UJ7Hq1auTmJhIaGgon376KY0bN2bPnj1069aN3bt3A3D//n2zBZbrxo0bAEZXstjb25Oens69e/ewtbV96nY2NjbExMQU+qU7IYR4rkroSGz06NFm60vxX20fHx+++OILdu7cyZgxY2jbti0xMTG8/37O5rG5z8vMLTcxPrnLMaB76TojI8MgiZnSzsLCQvH8qxBCFBcl9ZkY5MyARUZGEhcXR1paGhUqVKBZs2b4+voa/budF8VJbPDgwdy7d4+oqCicnJyoUqUKW7du5fTp0wDY2dkZPQ26sB5fyp8XYwnoadsJIURJUVKfid25c4f+/fvzxx9/UK1aNezt7bl8+TJ79uxh06ZNREREKD5nTHESu3v3Lh999JFuU+DSpUsTHh7O0aNHycjIoHnz5s9kJ4/cbfqzsrIMruWWGdvK/2nbCSFEiVFCk9iiRYu4ceMGGzZs0K2pAIiPj8ff35+QkBCmT5+uqC/FQ5Hu3bvzzjvv8PPPP+vKLC0teeONN+jSpcsz24qqWrVqANy8edPgWkpKCnZ2dgZHwxSmnRBClBSqbG2Bn+LowIEDjBkzRi+BATRv3hx/f3/279+vuC/FI7Hbt2+TkpJC+fLllUdqBnZ2djg5OXHmzBmDa2fPnqVBgwZmbSeEECVFSZ1OzMjIwMnJyeg1Jycn7ty5o7gvxSOxTp06AfDDDz8o7txcvLy8iI2N5cKFC7qyo0ePcunSJbp162b2dkIIURKosgv+FEd16tTh4MGDRq8dOHCAmjVrKu5L8UisTp062NjYMG/ePMLDw6lduzY2NjZ6y9JN2T7fFMOGDeO7775j4MCBDB48mKysLMLCwnBzc6NHjx5Azh6JnCJxigAAFv5JREFUp06domnTptSoUUNxOyGEKLHMPF2YnZ1NaGgokZGRpKWl4eHhwYwZM/JMKleuXGH+/PmcOHECjUZDo0aNmDhxIq+++mq+9xk8eDCffPIJGo2Gt956C3t7e27evMmuXbvYvn07M2fOVByz4iT2+L6JV69e1e3gkSt378RnkcQqVqxIeHg4QUFBLFmyBGtrazp16sSECRN0+x/Gx8czefJkgoKCdElMSTshhCipzD3SWrZsGREREcydO5cqVaqwcOFChgwZQnR0tN45kgDp6ekMHDiQOnXqsGbNGiwtLVm2bBn9+/dn165dVKpUKc/7dOvWjT///JOVK1eybds2ICeHWFlZ8eGHH/Lee+8pjllxEnNwcHhmG/wqUbt2bVatWpXndV9fX3x9fU1uZ0xew1whhChWtOYbianVatasWUNAQIBu+8Dg4GBat27Nnj17DHacP3z4MDdu3OC7777Tvac7f/58PD09OXDgQIGHJI8aNYp+/fpx+vRp0tLSKFeuHO7u7pQrV86kuBUnsSNHjpjUsRBCiGfLnCOxc+fOkZGRQcuWLXVlNjY21K9fnxMnThgksaZNm/LVV18ZbDSh1WoVL8yws7Mr9H67ipNYYGAgKpWK8ePHGwR94cIFFi5ciKOjI4GBgYUKSAghhDIqBSOxtLQ00tLSDMrt7Oz0XijOa6s+BwcHrl+/btDe0dERR0dHvbJ169aRlZVlNDH179+/wFhzqVQq1q1bp6iu4iQWGRmJSqVi9OjRBkns3r17HDx4EDs7O0liQgjxnKg0BSexdevWGT1iavTo0fj7++t+zszMBDBYL2BlZYVarS7wPnv27CEkJISBAwfi4uJicP3hw4cFPpL6448/SE9Px9LSssD75coziSUmJuLr66u3fZNWqzXY4V2r1erqyDZOQgjxHCl4JDZgwAB8fHwMyp/c1sna2hrIeTb2eCJTq9UFbgyxfv16goKC8Pb2znP7wYiIiDzbp6enM3fuXE6dOkWtWrWYO3duvvd7XJ5JrG7durz33nts2LAB+N9RLI8ePcqzszZt2ii+sRBCiEJSMJ345LRhXnKnBlNSUrCxsdGVp6SkULduXaNtsrOzmTNnDuHh4QwfPpxPPvnE5AWAx44dY8qUKSQnJ9O/f3/GjRtnsBIyP/lOJ3700UeULl0arVbL119/jUql4r333tPLyiqVCisrK1599VU6d+5sUvBCCCGenjm3lXJ1dcXGxoa4uDhq164N5IyQzp49S9++fY22mTlzJpGRkUyfPl13oolSDx48YP78+URERFC9enXWr19Ps2bNTI5bpdUqW6MZEBCASqUiMDBQ8e7CQgghnp3OrecUWGffj1MV9xccHMzmzZuZM2cOTk5OLFy4kCtXrrBz504sLS25ffs2tra2WFtbExMTg7+/PyNGjKBfv356/ZQpUybfDdZPnDjBlClTSEpKok+fPowfP96k41cepziJPe7OnTv88ccf3LlzBy8vLx48eKCbTxVCCPF8dG41u8A6+45OU9yfRqMhODiY7du3k5mZqduxo0aNGly9epWOHTsSFBSEr68vo0ePZt++fUb7GTFiBGPHjjUoV6vVLFq0iHXr1uHo6Mjnn3+ut6T/aZiUxBITEwkKCiI2Nla3Q8eZM2do3bo1/v7+9O7du1DB/BN0tuhV1CGYZF92ZImKuaTFC/8fs6XyHQiK2j7NNwC8WWFoEUei3PepYbxZcVhRh2GS72+btgmDMV6vf1ZgnZjY4rFi/JdffmHixIn8+eef9OrVi0mTJpnlJBHFS+wTExPp3bs39+/f5/G89+eff/L333/z6aefUqlSJXkuJoQQz4uCJfbFRZ8+fcjOzsbW1pYLFy4wbFje/9GhUqkIDw9X1K/iNfGLFi0iPT2dxo0bExISoisvV64cLVq0QKvVsnr1aqXdCSGEKCSVVlvgp7ho2rQpzZo1w8XFBQsLi3w/pqxwVDwSO3HiBCqVimnTplG5cmVdecWKFRk/fjzvvPMO58+fN+1bCSGEeHrZxfSsFSNyX9cyN8VJLPeNbWOP0LKysgB52VkIIZ6rkpPDnhnFWSf3JOQ5c+Zw+vRpXfnJkyeZM2cOKpUKN7f/a+/eY5o6/zCAP6XA8FbxgjoHbt6GFoaKWhSvMcxlFweSuGWJeEVNpsExEwRnzP5wc8YYN38Y3OZUdIuLbhZnJMFkXhbFeGNkKu6mbGKmgiIqoFjb8/uDtFChtEXoOd/yfJJmoT3n7dPE7Jv3ct43qu0TEhFRs3Q2m9uXv/O4J5aWloYFCxaguLgYxcXFjjHL2bNnQ1EU6PV6LFmypN2CEhHRUzQ056UWj3tiJpMJ2dnZ6NevHxRFcXoNGDAAX3zxBcaPH+91gLKyMixbtgwmkwkmkwkZGRmorKxss/ta0/7q1auRkpLi9W8hIvIpRXH/8nMe98QAYOrUqZg8eTIuXryIa9euQa/XIzw8HFFRUa2aD7t79y7mzp2Lx48fIzU1FVarFd988w3++OMP7Nu3z+Xpy57e15r29+3bh3379sFkMnn9e4iIfMmTXez9ncsilpWVBZ1Oh6ysLKejVwICAhATE4OYmJhn/vKdO3fi5s2bOHjwIAYPHgwAGDFiBObPn4+8vDyXJ4N6ep837VutVuTk5DR7ZAERkSZ1gJ6WOy67T2azGWazGbW1te325YcOHYLJZHIUGACIj4/HwIEDcejQoWe+z9Pr6urqMHPmTPzvf/9DYmJik0PhiIg0yWpz//Jzqq2Jv3fvHsrKyppd0RgVFYWLFy8+033etF9XV4fq6mps2rQJ69evR2CgV6OsRETq4JyYd3NibcnVUdgAEBYWhurqajx48KDJKdKe3udN+127dsXhw4dZvIhIlg7Q03LH7f+1ly9f7nKBxdN0Oh1yc3M9urampgYAmt1+334gWm1tbZMi5ul93rRv3+qEiEgUhUXMbRErLi5224hOp3Psau8pmwcP4TVXWDy9r7XtExGJwZ6Y+yJmNBoRFBTU5l9sPzDNvmVVY/b3mjtUzdP7Wts+EZEYHWDOyx23RSwnJ6ddVuv1798fAFBRUdHks/LychgMhmbPmvH0vta2T0QkRgfYVsod1VYyGAwGhIeH49KlS00+KykpcezV2Nr7Wts+EZEYLGLqLbEHgOnTp+PUqVO4cuWK473CwkKUlpbijTfeeOb7Wts+EZEIXGLvuieWlJQEnU7XrkNuixYtwoEDBzBv3jwsWLAAdXV12LZtG6KiopCYmAigfu/DoqIixMbGIiIiwuP7vLmOiEgkLuxwXcQ+++yzdv/ynj174ttvv8W6deuwefNmhISEICEhARkZGY5l/WfPnkVWVhbWrVvnKGKe3OfNdUREEilcYq/enJjdoEGD8PXXX7v8PDk5GcnJyV7f5+11jR05csSr64mIVMGemPpFjIiIWokLO1jEiIikUqxWtSOojkWMiEiqDrD60B0WMSIiqdgTYxEjIpJKsbEnxiJGRCQU58QAnaJwUJWIiGTiWSRERCQWixgREYnFIkZERGKxiBERkVgsYkREJBaLGBERicUiRkREYrGIERGRWCxiREQkFosYERGJxSJGRERisYgREZFY3MWeSJjq6mo8fPgQAQEB6NKlC0JCQtSO1CJpeQGZmTsqFjEijXvy5AnMZjPy8/Nx8eJFVFdXO31uMBgQHR2NN998EzNmzEBQUJBKSetJywvIzEz1eBSLBsyZM8fre3Q6HXJzc9shjXvS8gIyMwNAZWUlFi5ciMuXL2PgwIEwGo3o06cPnnvuOQBAXV0dysvLUVJSgtLSUrzyyiv46quv0KNHD+b148zUgD0xDRgwYAB++OEH6HQ69OnTB3q9Xu1ILZKWF5CZGQDWr1+P69evY8eOHRg/fnyL1xYWFmL58uXYsGEDPv30Ux8ldCYtLyAzMzWikCbs3btXGTZsmJKZmal2FI9Iy6soMjObTCZl69atHl+fk5OjxMfHt2OilknLqygyM1MDrk7UiFmzZmHp0qXIy8tDQUGB2nHckpYXkJkZAAIDPR8wCQkJwcOHD9sxjXvS8gIyM1M9FjENef/99zF06FBs2LABNptN7ThuScsLyMs8duxY5ObmorS01O21paWl2L59O8aOHeuDZM2TlheQmZkacGGHxjx+/Bh1dXXo1q2b2lE8Ii0vICvz9evXMWfOHJSXl8NkMiEqKgp9+/ZFSEgIdDodHj16hIqKCpSUlKCwsBDdu3fHd999h5deeol5/TgzNWARI9K4qqoq7NixA/n5+SgrK2v2mvDwcEyfPh2pqano2bOnjxM6k5YXkJmZ6rGIaZS0hy2l5ZXq/v37uHXrFmpqamCz2dClSxf0799fs71KaXkBmZk7MhYxjZD2sKW0vN4qKyvD+fPnkZSUpHYUImoBi5gGSHvYUlre1vjpp5+wcuVKXL58We0oXpFWfKXlBWRm9mtqre2nBhkZGcqYMWOUwsJCt9eePHlSGTNmjJKVleWDZM2Tlrc1rl27puzfv1/tGF47cOCAMmzYMLVjeExaXkWRmdmfsSemAXFxcViwYAGWLFni0fVbt27F7t27cfLkyXZO1jxpeV3xx3m8srIynDt3DjNnzlQ7ikek5QVkZvZn3HZKI6Q9bCktL+D/83gAEBERgYiICLVjeExaXkBmZn/GIqYB9octp02bhoEDB7Z4rRYetpSWF2g6jzd58mSX83irVq3Cnj17xM3jaZXEHq/EzB0VhxM1QNrDltLyAsDKlStx5MgRbN682eNNXl999VVu8toKEnu8EjNTPRYxjZD2sKW0vFLn8aQdISNx5arEzNSAw4kaERoaivT0dKSnp4t42FJaXkDmPJ60I2QkHmsiMTM1oubSSHLtwYMHSnl5uXL79m3l4cOHascRb+nSpcqkSZOUq1evur326tWryqRJk5TFixf7IJl7ko6QkXisicTM1IA9MY3w5zF5LTwcmpmZiTlz5mDGjBkez+NlZWWplrexWbNm4datW9iyZQumTp2K1157Te1ILZLY45WYmepxTkwD/H1MXiu7X0ibx2vMZrMhKSkJtbW1OHz4MAICtHmK0rJly/Dbb78hNzfXo5Wrc+fOxfDhw/Hll1/6KGFTEjNTAxYxDfD3lXNafDhUyjxeYxKOkJG4clViZmrAIqYBUlfONcbnashOYo9XYmaqxzkxjZA2Ju/Pc3j0bCSuXJWYmeqxJ6YB0sbk/X0Oj4gjC3KwiGmAtDF5f5/Do46HIwtysYhphKQxeX+Yw6P2IW2HEYAjC9JxTkwjpI3JS5vDI9+QtsMIwB07pGMR0yCDwQCDwaB2DJck7mJPvrF27VqMGDECa9asQXx8PNatW6d2JLeOHTuG1NRUtwUMAOLj47Fw4ULs3r3bB8nIEyxiAqm9A4bk3S+o/UnbYQTgyIJknBMTSAs7YEiawyPfk7LDCCBvdTA5YxETSGs7YEiYwyPfk7DDCCBvdTA5YxHTKD6nQuQ7HFmQi0VMI/icCpE2cGRBFhYxDeBzKkRErcPViRrA51SIiFqHRUwD+JwKkXok7jJCDVjENILPqRCpQ+IuI9SARUwDuAMGkXok7jJCDbiwQwP4nAqR+rKzs7FlyxZ8/vnnInYZoXosYhrB51SI1CVplxFqwCKmQXxOhUgdUnYZoQYsYkREJBb7y0REJBaLGBERicUiRqKlpKQgMjIS77zzTpP3Gr+io6MRFxeHlJQUnD171nFtZmZmk2ube2VmZjp9r8Viwe7du/Huu+8iLi4O0dHRmDhxItLS0nDu3DmvfsOWLVsQGRmJ48ePe/37n85/8OBBp8/PnTvn9HlKSgoA4OrVq4iMjMRHH33k9XcSaQmLGHUIFosFVVVVOHPmDObPn+91oWns1q1beO+997B27VoUFxejqqoKFosFFRUVKCgowOzZs5GXl+dRW+Xl5di2bRv69u2LiRMntjqTXeMCDcDl7xw0aBBiY2Oxf/9+Vc+lI3pWLGLkt4xGI44fP46jR4/i0KFDSEtLg06ng8ViwaZNmwAAWVlZOH78uONlNBoBACNHjnR6334ytdVqxQcffIALFy5Ar9cjNTUV+/fvR35+PlatWoXOnTtDURSsXr0a//zzj9uMO3fuRG1tLZKSktpkp4gzZ844/f10UWssOTkZNpsNOTk5z/y9RGrhjh3kt4KCgtCvXz/H30OGDMGFCxdw9OhRFBUVwWKxoHv37ujevbvTPQAQHBzsdK9dXl4eioqKAAAff/yx0zDm4MGD0atXL6xYsQIWiwVmsxnp6eku8z169Ag//vgjACAhIeGZfuvzzz+PGzduoLS0FBUVFQgLC4PVasWvv/4KAOjfvz/+++8/p3umTZsGnU6Hn3/+2XEPkTTsiVGHMnToUAD1D7bevXvX6/vtw4QDBgzArFmzmnz++uuvY8WKFdi+fTsWLVrUYlunT59GVVUVevbsiejoaMf79nmuCRMmOF3v6n0AiIiIcBQhe++rpKQENTU1CAoKQkxMTJN7evXqBaPRiCdPniA/P9/NLyfSJhYx6lD+/PNPAIBer0doaKjX9xcXFwOo3+9Sp9M1+Vyv12Px4sWYMGECunbt2mJbhYWFAICoqKg22R0iNjYWQEMRs8+HDR8+HJ06dWr2HnvxfHoYkkgKDieS37JYLLh58yYURcH9+/dRUFCAY8eOAQDi4uIQHBzsVXuVlZV4/PgxADQZemv8mZ1er29xiM4+1DdkyBCvcrgSGxuLgoICRxGz/zc2Nhb37t1r9h57z9SehUgaFjHyWyUlJZgyZUqT9zt37oyMjAyv22vcW7JarU6fffjhhzh16pTTey+88AKOHDnisr3bt28DaFoQW8veE/v7779RWVmJ8+fPAwBGjx7tMkfv3r0BAHfu3IHFYnHMCRJJweFE6hCCg4PRt29fJCQk4Pvvv8fw4cO9biM0NBQhISEAgBs3bjxzpsrKSgBAly5dPLre3Q5xRqMRnTp1gqIo2LNnD6qqqgA0FLfmNP5u+/VEkrAnRn5rxIgR2Lt3b5u2OWbMGJw4cQKnT5+G1Wp1LIvfuXOn45qUlBSP5pjsRcnV0nqbzeb099PDlU8LDAxETEwMTp8+jV27dgEAXnzxRUdvqzmN5/W8OZiVSCvYEyPywowZMwAAFRUV2L59e5PPrVYrampqPGqrR48eANDkevuQXm1trdP7N2/edNumvdflSS8MAB48eACgvpg1ftSASAoWMSIvJCYmYtSoUQCAjRs34pNPPsGlS5fw77//Ij8/H8nJybh06ZJHbYWHhwNoWpzsc2T258isVivMZrNHiy+eLlqjR49u8fo7d+4AqO+x8fwskoj/aom8oNPpkJ2djVGjRkFRFOzatQvJycmYPn060tPT8fvvvwMADAYD0tLSWmxr5MiRAIArV644vZ+QkOAoKKtWrUJUVBQyMzMdu4m0ZNSoUU7FyF1P7K+//gJQv8yfSCIWMSIv9e7dG7t27cL69esxbtw49OnTB0FBQQgLC8OkSZOwZs0aHD16FElJSS22YzKZAABFRUVO819GoxEbN27E4MGDERwcjMjISGzatAlvv/2222zdunVzLNkPDQ3FoEGDWrzeXnTHjx/vtm0iLeKhmEQqefLkCaZMmYLbt2/DbDZ71NNqS9XV1Rg3bhx0Oh1OnjwJg8Hg0+8nagvsiRGpJDAwEImJiQCAEydO+Pz7f/nlF1gsFkybNo0FjMRiESNS0fz589G5c2ePj25pS2azGQEBAVi6dKnPv5uorbCIEakoLCwM8+bNw5UrV3y6f2FZWRlOnDiBt956Cy+//LLPvpeorXFOjIiIxGJPjIiIxGIRIyIisVjEiIhILBYxIiISi0WMiIjEYhEjIiKx/g+n+DMbWWeuLAAAAABJRU5ErkJggg== " >
In [48]:
x0 = {"dna_pttr_hrpR":1,
      "dna_constit_ttrS":1,
      "dna_constit_ttrR":1,
      "dna_constit_lacI":1,
      "dna_plac_hrpS":1,
      "dna_hrpL_gfp":1,
      "dna_constit_mScarlet":5,
      "phosphate_P":100,
      "protein_RNAP":15,
      "protein_RNAase":30,
      "protein_Ribo":120,
     "protein_lacI": 0}
In [49]:
IPTG_um_list = np.logspace(-3, 2, 6).round(decimals=6)
tt_um_list = np.logspace(-3, 2, 6).round(decimals=6)

fig, heat_map = simulate_heatmap('ligand_tt',tt_um_list,'inducer_IPTG',IPTG_um_list,x0=x0,normalize=True,title='Compiled Full AND gate \n No lacI Spike', CRN=CRN_extract_1)

heat_map.set(xlabel='IPTG (uM)', ylabel='Tetrathionate (uM)')
heat_map.invert_yaxis()
plt.yticks(fontsize=18, fontstyle='normal', rotation = 0)
fig.savefig("./figures/fig6/fig6b.png")
plt.show()

Early Circuit Dynamics in no lacI Spike Model

In [50]:
timepoints = np.linspace(0, 200, 1000)

x0 = {"dna_pttr_hrpR":5,
      "dna_constit_ttrS":1,
      "dna_constit_ttrR":1,
      # Cole origin
      "dna_constit_lacI":40,
      "dna_plac_hrpS":1,
      "dna_hrpL_gfp":1,
      "dna_constit_mScarlet":5,
      "phosphate_P":100,
      "protein_RNAP":10,
      "protein_Ribo":90.,
      "ligand_tt": 1,
      "inducer_IPTG": 0.,
      "protein_lacI": 0,}
Re1 = CRN_extract_1.simulate_with_bioscrape(timepoints, initial_condition_dict = x0)

results_melted = pd.melt(Re1, var_name='species', id_vars=['time'])
results_melted.head()
Out[50]:
time species value
0 0.000000 complex_2x_phosphate_P_protein_ttrR_dna_pttr_h... 0.000000e+00
1 0.200200 complex_2x_phosphate_P_protein_ttrR_dna_pttr_h... 5.330060e-13
2 0.400400 complex_2x_phosphate_P_protein_ttrR_dna_pttr_h... 2.037818e-11
3 0.600601 complex_2x_phosphate_P_protein_ttrR_dna_pttr_h... 1.643681e-10
4 0.800801 complex_2x_phosphate_P_protein_ttrR_dna_pttr_h... 7.002289e-10
In [51]:
Re1 = CRN_extract_1.simulate_with_bioscrape(timepoints, initial_condition_dict = x0)
/Users/lianamerk/anaconda3/envs/python3.7-bioscrape/lib/python3.7/site-packages/biocrnpyler-0.1-py3.7.egg/biocrnpyler/chemical_reaction_network.py:870: UserWarning: The following species are uninitialized and their value has been defaulted to 0: inducer_IPTG, 
In [52]:
p = bokeh.plotting.figure(title = 'Early Circuit Dynamics, pLac', x_axis_label = 'Time (min)', y_axis_label = 'Simulated Concentration (uM)', height = 2000, width = 4300)

p.line(timepoints, Re1['complex_dna_plac_hrpS_protein_RNAP'], legend_label = 'Active pLac', color = 'Green', line_width = 14)
p.line(timepoints, Re1['rna_constit_lacI'], legend_label = 'LacI RNA', color = 'orange', line_width = 14)

p.line(timepoints, Re1['protein_hrpS'], legend_label = 'hrpS protein', color='lightsteelblue', line_width = 14)

p.line(timepoints, Re1['complex_dna_hrpL_gfp_protein_hrpR_protein_hrpS_protein_RNAP'], legend_label = 'hrpS, hrpR both bound to pHRPL', color='darkorchid', line_width = 14)
p.line(timepoints, Re1['rna_hrpL_gfp'], legend_label = 'GFP RNA', color='navy', line_width = 14)

p.add_layout(p.legend[0], 'right')
p.legend.title = 'Species'


p.legend.glyph_width = 100
p.legend.padding = 100
p.legend.title_text_font_size = '80pt'
p.legend.label_standoff = 50
p.legend.glyph_height = 100

p.xaxis.major_tick_line_width = 10
p.yaxis.major_tick_line_width = 10

p.xaxis.minor_tick_line_width = 10
p.yaxis.minor_tick_line_width = 10
p.axis.minor_tick_out = 20
p.axis.minor_tick_out = 20

p.axis.major_tick_out = 40
p.axis.major_tick_out = 40
p.yaxis.axis_label_standoff = 50
p.xaxis.axis_label_standoff = 50

p.xaxis.axis_line_width = 8
p.yaxis.axis_line_width = 8

p.yaxis.major_label_standoff = 30
p.xaxis.major_label_standoff = 30

p.ygrid.grid_line_width = 10
p.xgrid.grid_line_width = 10

p.axis.axis_label_text_font_size = "80pt"
p.legend.label_text_font_size = "80pt"
p.title.text_font_size = '100pt'
p.yaxis.major_label_text_font_size = "80pt"
p.xaxis.major_label_text_font_size = "80pt"

#p.y_range = Range1d(0, 0.9)
bokeh.io.show(p)

Computing Environment

In [53]:
%load_ext watermark
%watermark -v -p numpy,pandas,scipy,bioscrape,biocrnpyler,bokeh,jupyterlab,holoviews
CPython 3.7.6
IPython 7.13.0

numpy 1.18.1
pandas 1.0.3
scipy 1.4.1
bioscrape 1.0.0
biocrnpyler unknown
bokeh 2.0.1
jupyterlab 2.1.0
holoviews 1.13.2